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ABSTRACT 


Lieu et al. (2015) have recently claimed that it is possible to substantially improve the sensitivity of radio-astronomical 
observations. In ess ence, their propo sal is to make use of the intensity of the photon shot noise as a measure of the 
photon arrival rate. Lieu et al. ((2015) provide a detailed quantum-mechanical calculation of a proposed measurement 
scheme that uses two detectors and conclude that this scheme avoids the sensitivity degradation that is associated with 
photon bunching. If correct, this result could have a profound impact on radio astronomy. Here I present a detailed 
analysis of the sensitivity attainable using shot-noise measurement schemes that use either one or two detectors, and 
demonstrate that neither scheme can avoid the photon bunching penalty. I perform both semiclassical and fully 
quantum calculations of the sensitivity, obtaining consistent results, and provide a formal proof of the equivalence of 
these two approaches. These direct calculations are furthermore shown to be consistent with an indirect argument 
based on a correlation method that establishes an independent limit to the sensitivity of shot-noise measurement 
schemes. Furthermore, these calculations are directly applicable to the regime of interest identified by Lieu et al. 
Collectively, these results conclusively demonstrate that the photon-bunching sensitivity penalty applies to shot-noise 
measu rement schemes just as it does to ordinary photon counting, in contradiction to the fundamental claim made by 
Lieu et al. (20151. The source of this contradiction is traced to a logical fallacy in their argument. 


Subject headings: instrumentation: miscellaneous 


1. INTRODUCTION 

In the infrared, optical, or x-ray bands, detection sensi¬ 
tivities are ultimately limited by the Poisson statistics of 
photon counting, with r.m.s. count fluctuations given by 
\/N where N is the mean number of photons collected 
Gehrels 1986). Thus the Poisson uncertainty in the 
ux measured tor an astronomical source is proportional 
to the square root of the total intensity of the radia¬ 
tion falling on the detector. Meanwhile, sensitivities for 
radio-astronomical o bservations a re calculated using the 
radiometer equation (Dicke 1946), which states that the 
measurement uncertainty is proportional to the the total 
radiation intensity rather than its square root. 

A transition between these two regimes - radio and op¬ 
tical - is therefore inevitable, and corresponds to a shift 
from a classical description involving fields and waves 
to a quantum description involving photons, som eti mes 
referr ed to as the radio-optical dichotomy (Nityananda 
19941 and ultimately stemming from the wave-particle 
duality of quantum mechanics. The nature of this tran¬ 
sition was clarified through the demonstration of corre¬ 
lated photon arrivals at two independent detectors by 
Hanbury Brown & Twiss (1956; HBT), using a setup 
similar to that illustrated in Figure |T] The HBT cor¬ 
relations are a manifestation of photon bunching, which 
causes the photon arrivals to be clustered in time rather 
than being purely random. As described in more detail 
in section [2j bunching causes the photon count fluctu¬ 
ations for a single detector to increase to y/N(\ + n) 
rather than the usual y/N for Poisson statistics. Here 
n represents the photon mo de occupation number for a 
detector with unit efficiency (Zmui dzinas|2003 ), given by 
the Bose-Einstein formula n = 1 j(^ hv l kl — 1) for thermal 
blackbody radiation at a temperature T. Bunching is 
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Fig. 1. — A two-detector exp eriment similar to that used by 
Hanbury Brown & Twiss| fl956| and others (e.g., Harwit 1960) 
to demonstrate photon correlations. A bright thermal light source 
with narrow spectral bandwidth A v illuminates two photon detec¬ 
tors via a 50/50 beamsplitter, producing photocurrents Ii(t) and 
I 2 (t) • Photon bunching results in a nonzero correlation between 
the photocurrents: (Ii(t)l 2 (t)) ^ 0. The unused input port of the 
beamsplitter is terminated with a cold (dark) absorber to prevent 
stray light from entering. 


usually ignorable for astronomical observations made in 
the infrared to x-ray bands because n « lQin contrast, 
photon bunching is a large effect in the radio band since 
n ^ kT/hu » 1. Furthermore, at radio wavelengths 
both N and n scale with the intensity of the radiation be¬ 


ing detected; therefore, y/N( 1 + n) is also proportional 
to the intensity, in agreement with t he Dicke equation. 
In a recent paper, Lieu et al. (2015) claim to have found 


1 Thermal radiation at optical wavelengths with high occupa¬ 
tion number n may readily be generated in the laboratory using 
stochastically modulated coherent laser radiation, e.g. p roduced 
by scattering from rotating ground glass (R ousseau 1971). 
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Fig. 2.— A simple setup consisting of a thermal radiation source, 
an optical bandpass filter with transmission bandwidth Au, and an 
ideal photon detector. The detector output is represented by the 
photocurrent I(t). The thermal source and filter comprise the light 
source in Figure [T] 



Fig. 3.— The photocurrent noise spectrum Sj(u) for an ideal 
photon detector illuminated with filtered thermal radiation consists 
of three components: (1) a DC term contributed by the mean 
photocurrent; (2) a photon bunching component that extends to 
±Av, where Au is the bandwidth of the radiation; (3) a white shot 
noise term that rolls off at ± u d , the detector bandwidth. Both the 
DC photocurrent and the shot noise intensity are proportional to 
the average photon arrival rate F = nAu. The shaded region shows 
the portion of the spectrum that is available for measurement of the 
shot noise intensity without interferenc e fro m the bunching noise 
component. See section [ 2 ] and equation (|20|) for details. 


a method for avoiding the extra noise associated with 
photon bunching and thereby potentially increasing the 
sensitivity of radio telescopes by \J 1 + n, which is a large 
factor, e.g., over an order of magnitude for the example 
of Arecibo described in their paper. Such a possibility is 
of obvious interest given the large sums spent on the con¬ 
struction of radio telescopes and the associat ed receiving 
equipment. The essence of the Lieu et al. (2015) pro¬ 
posal is to use the wide-bancl shot noise at the output 
of a fast photon-counting detector as a measure of the 
radiation intensity^ Consider the simple single-detector 
setup illustrated in Figure [2] It is helpful to visualize 
the noise spectrum (the power spectral density, or PSD) 
at the output of the detector as illustrated in Figure [3] 
which graphically summarizes the quantum-mechanical 
calculatio ns presented later in section [2] The shot-noise 
spectrum (Schottky 1918) is white andieatureless within 
the output bandwidth of the detector, and has an inten¬ 
sity that is proportional to the mean photon arrival rate 


2 A fast photon detector operating at radio frequencies may rep¬ 
resent a serious technical challenge, but not one of fundamental 
principle: tunnel junction detectors offer one possible method of 
implementation (Tucker &; Millea 1978| |SchoeIkopf et al.|1999|. 


T. Meanwhile, the bunching noise component is confined 
to lower frequencies, determined by the bandwidth Av 
of the radiation being detected. In principle, the radia¬ 
tion bandwidth Av may be made arbitrarily small using 
narrow-band filters preceding the detector, so the use of a 
fast detector with an output bandwidth v ( j >> Av allows 
the region of the photocurrent noise spectrum where the 
white shot noise dominates to be accessed and measured 
with appropriate signal processing techniques. Clearly, 
it is advantageous to use a large measurement bandwidth 
B = v ( j ~ Av, since the fractional precision with which 
the shot noise intensity may be determined cannot be 
better than 1 / \/BT, where T is the measurement time 
(Di cke||1946D . 



Fig. 4.— The noise spectrum Sj A (V) for the output 7 a of the 
two-detector scheme proposed by|Lieu et al. | (|2015|) and illustrated 
in Fig. Here 7 a = 7i — I 2 represents the difference in the pho- 
tocurrents for the two detectors. Taking this difference eliminates 
the DC component as well as the bunching component in the spec¬ 
trum, leaving only the white shot^ noise component that is propor¬ 
tional to the mean photon rate F. The full spectrum is available 
for measurement of the shot noise intensit y, as i llus trated by the 
shaded region. See section |3| and equations < |39[ ), ( |40[ ), and ( |41| > for 
details. 


Alternatively, as specifically proposed by |Lieu et ak| 
(20151 and shown in Figure [l] a 50/50 beamsplitter or 
its radio equivalent may be used to feed two detectors. 
Each detector individually has an output noise spectrum 
similar to that shown in Figure [3] although with half 
the total photon rate (T/2) per detector. The DC term 
may be eliminated by taking the difference of the two 
photocurrents. Differencing also eliminates the bunching 
noise lying in the frequency interval [— Av, Av], because 
this component is fully correlated at the two detectors, as 
is demonstrated through a quantum-mechanical calcula¬ 
tion in section [3] Indeed, this component is responsible 
for the HBT correlations. Thus, only the white shot- 
noise spectrum survives after taking the difference, as 
shown in Figure HI the shot noise intensity may then be 
measured using relatively simple signal processing tech¬ 
niques. Although the measurement bandwidth may now 
be increased to B = v d instead of v d — Av, the res ult¬ 
ing im provement is modest when v d » Av. Lieu et al. 
(20151 present a full quantum-mechanical calculation of 
the sensitivity of such a shot noise measurement scheme 
using two detectors, which is a nontrivial task involving 
computation of eighth-order moments of photon oper¬ 
ators, and conclude that the \fN Poisson uncertainty 
may be achieved instead of the usual bunching-degraded 


vW + n) uncertainty as expressed by the Dicke equa¬ 
tion. This result is quite surprising, and if correct and 
amenable to practical implementation, would represent 
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a significant discovery with the potential to stimulate 
large advances in rad io astronomy. However, as demon¬ 
strated in section m the fundamental conclusion of the 
Lieu et al. (2015) paper rests on a logical fallacy and 
is therefore not valid. Section [TO] also contains a simple 
intuitive argument that demonstrates why the measure¬ 
ment scheme proposed by Lieu et al. is in fact subject to 
the photon bunching penalty; those uninterested in the 
detailed calculations in the following sections may wish 
to jump straig ht to section [lOl 
The work of Lieu et al. (2(jl5|, and the quantum cal¬ 
culations presented here, may pose a challenge to those 
who are more familiar with classical concepts such as 
fields and voltages than with photon operators and quan¬ 
tum mechanics. Nonetheless, the essence of the problem 
is quite easy to understand by use of a familiar anal¬ 
ogy. The analogy relies on the fact that thermal photon 
bunching can be correctly described by a photon arrival 
rate that varies with time in a random way, as will be 
discussed below. Imagine listening to the sound of rain 
landing on a roof: this is the shot noise produced by 
the random arrivals of a large number of individual rain¬ 
drops. The intensity of the sound depends on how hard 
it is raining, i.e. the arrival rate of the raindrops]^] If 
the raindrop arrival rate changes with time, as often oc¬ 
curs over timescales of seconds to minutes, the inten¬ 
sity of the sound will vary accordingly. Thus, while a 
measurement of the total precipitation may be made by 
integrating the acoustic shot noise intensity over time, 
this shot noise measurement will reflect the fluctuations 
of the raindrop arrival rate in the same way as would 
a direct measurement, e.g. observations of the water 
level in a standard rain gauge. Note that the spectral 
character or the “sound” of the acoustic rain noise re¬ 
mains constant as the intensity changes. Furthermore, 
note that the connection between the two measurement 
methods - acoustic noise vs. rain gauge - is purely classi¬ 
cal and has nothing to do with quantum mechanics. This 
statement is also true for the photon detection problem. 
Indeed, the output of a photodetector is an entirely clas¬ 
sical quantity - a train of electrical pulses - whose prop¬ 
erties are fully specified by the statistics of the photon 
arrival times. While we may need to turn to quantum 
mechanics to calculate the arrival time statistics, once 
the arrival time statistics are known, in principle we can 
generate a simulated classical pulse train numerically, as 
illustrated in Figure [5] and use this time stream to cal¬ 
culate any other quantity of interest, e.g. the mean and 
variance of the photon counts, or the mean and vari¬ 
ance of the shot noise intensity, or the correlation be¬ 
tween the photon counts and shot noise intensity, etc. 
These quantities are all related to various moments of 
the same classical time stream. It is therefore not sur¬ 
prising that bunching affects standard photon counting 
measurements and photon shot noise measurements in 
similar ways, and therefore the sensitivity degradation 
due to bunching cannot be avoided. Indeed, in section [8] 
I present a calculation that demonstrates that the shot 
noise intensity has a nonzero correlation with the photon 
counts, and then use this correlation to establish a rig¬ 
orous sensitivity bound for the shot noise measurement. 
This bound shows that the shot noise measurement is 


5 The raindrop size should be kept fixed for the analogy to hold. 


subject to the same \/l + n sensitivity degradation due 
to photon bunching as for ordinary photon counting. 



Fig. 5. — Top left: Simulated output time stream for a single de¬ 
tector (see diagram in Figure[2]l illuminated by a coherent source or 
a thermal source with low occupation number. The photon arrival 
rate is constant with time, which is the case of Poisson statistics. 
Bottom left: Differenced output of a beamsplitter-fed pair of de¬ 
tectors (Figure Q for the same Poisson case. Top right: Output 
time stream for a single detector illuminated with strongly bunched 
thermal radiation (high occupation number). Bottom right: Dif¬ 
ferenced output of a detector pair for the bunched case. Horizontal 
and vertical scales are in arbitrary units. 


Let us continue to accept the claim that photon bunch¬ 
ing can be correctly described by a photon arrival rate 
that varies randomly with time. Would we expect to see 
the noise spectra illustrated in Figures [3] and |dj which 
were derived from the quantum calculations presented in 
sections [2] and [3] ? It is helpful to visualize the detec¬ 
tor outputs as a function of time as shown in Figure [5j 
The detector bandwidth is assumed to be larger than 
the photon arrival rate, u,i > T, so the detected photons 
are visible as sharp, well-separated output pulses. For 
the case of two detectors, taking the difference of the 
two outputs means that the pulses may be positive or 
negative depending on which detector receives the pho¬ 
ton. The two subplots on the left correspond to the case 
that the photon arrival rate is kept constant, which pro¬ 
duces Poisson statistics; in contrast, the two subplots on 
the right were generated using a time-variable photon 
arrival rate in order to simulate photon bunching. Imag¬ 
ine that these output time streams are averaged over 
a timescale r that is long compared to 1/T. For the 
single detector case, Poisson arrivals would produce a 
DC component along with small fractional fluctuations 
of order 1/vTv. Meanwhile, a time-variable arrival rate 
would lead to a DC component along with significantly 
larger fluctuations, in accordance with the noise spec¬ 
trum shown in Figure [3] For two detectors, the posi¬ 
tive and negative pulseslargely cancel when performing 
the time average. This cancellation occurs for both the 
Poisson and bunched cases, in agreement with the noise 
spectrum shown in Figure[4] However, both positive and 
negative pulses deliver the same high-frequency energy, 
on average, to the subsequent circuitry and thus con¬ 
tribute equally to the shot noise intensity. Therefore, 
the shot noise intensity for the single-detector and two- 
detector cases should be the same. Thus, our conclusion 
is that a description of photon bunching in which the 
photon arrival rate varies randomly with time could in¬ 
deed reproduce the noise spectra in Figures [3] and |4j 
Is it in fact correct to view photon bunching as result¬ 
ing from a time-varying photon arrival rate? Is Figure [5] 
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a faithful depiction of t h e photo n bunch ing? Indeed, 
this was how Hanbury Brown & Twiss (1957) viewed 
the phenomenon in their original work. In their words: 

we are dealing essentially with an interference phe¬ 
nomenon which can be interpreted, on the classical wave 
picture, as a correlation between intensity fluctuations 
due to beats between waves of different frequency; the 
concept of a photon need only be introduced at the stage 
where energy is extracted from the light beam in the 
process of photoemission.” This physical picture is why 
the excess noise due to photon bunching continues to be 
referred to as “wave noise”. |Hanbury Brown fe Twiss| 


(1957) computed the effect using exactly this semiclas¬ 
sical approach, in which the light is first treated as a 
classical wave, consisting of a random superposition of 
components at different frequencies, resulting in an in¬ 
tensity that has fractional variations of order unity that 
occur on a “coherence” timescale r ~ Az/ -1 that is set 
by the fastest beat frequency that can be produced if the 
spectrum is restricted to an optical bandwidth Aza The 
photoemission rate is assumed to be proportional to the 
light intensity, and therefore the output of each photon 
detector may be described by a compound Poisson pro¬ 
cess in which the photon arrival rate varies stochastically 
with time. The classical light intensities calculated for 
the two detectors sho wn in Figure [T] would be identical; 
Hanbury Brown & Twiss (1957) therefore conclude that 
although the photoemission rates for the two detectors 
both fluctuate, the fluctuations of these rates are per¬ 
fectly correlated, and this leads to a nonzero correlation 
of the detector outputs. A similar semiclassical approach 
involvin g a co mpound Poisson process was described by 
Mandel 1959]). Section [5| presents a semiclassical analy¬ 
sis ol the sensitivity of a shot-noise measurement scheme 
using a single detector; the case of multiple detectors in 
treated in section El The conclusion of the semiclassi¬ 
cal analysis for both cases is that the shot noise schemes 
cannot improve on the y/N(l + n) bunching-limited sen¬ 
sitivity for standard photon counting. 

In addition to this historical basis, the interpretation of 
photon bunching in terms of a time-variable photon ar¬ 
rival rate is both supported by experiment and fully con¬ 
sistent with the predictions of quantum mechanics. The 
fact that the photon arrival rate for thermal radiation 
does indeed vary wi th time was directly demo nstrated 
in the laboratory by Morgan & Mandel (19661 through 
measurement of the correlation of the arrival times of in¬ 
dividual photo ns at a single detector. Indeed, Morgan & 
Mandel| ( 1966|) give a simple, concise description of the 
photon bunching effect: “In time intervals of order or less 
than the coherence time of the light, the probability of 
counting two pulses is greater than that expected for ran¬ 
dom events”, just as depicted in Figure [5} Furthermore, 
the results of the semiclassical analysis which invokes a 
stochastic, time-variable photon arrival rate are repro¬ 
duced by a full quantum calculation. |Kelley fc Kleiner| 
(1964) describe a quantum-mechanical theory of photon 
detection that uses a density matrix to describe the state 
of the electromagnetic field; the HBT photon bunching 
effect may be studied through use of a density matrix ap¬ 
propriate for thermal radiation. A fully quantum analy¬ 
sis for shot-noise measurements using a single detector is 
described in section [7J and agrees with the semiclassical 


results in section [5j As discussed in section [TJ the ex¬ 
tension of the fully quantum analysis to the case of two 
detectors is straightforward and agrees with the corre¬ 
sponding semiclassical analysis in section [6] Thus, both 
the semiclassical and fully quantum calculations show 
that the bunching noise cannot be evaded through use 
of a shot noise measurement scheme, whether one uses 
one or two detectors. This agreement is a reflection of 
the equivalence of the quantum an d sem iclassical descrip¬ 
tions of light as shown by Sudarshan (1963), who made 
use of the coh erent state representation introduced by 
Glauber (19631. In Appendix |F| the equivalence of the 
quantum and semiclassical (i.e., compound Poisson) de¬ 
scriptions of the photocurrent statistics is demonstrated 
explicitly. Thus, the interpretation of photon bunching 
in terms of a time-varying photon arrival rate as illus¬ 
trated in Figure[5]is in fact predicted by the full quantum 
theory and also supported by laboratory measurements. 

A potentially confusing aspect of the discussion is the 
fact that the shot noise spectrum is white - indeed, the 
output spectrum for the two-detector case (Figure 0]) is 
flat and featureless. Where is the bunching noise hiding? 
The answer is simple: the shot noise spectrum is white 
regardless of whether the photon arrival rate is constant 
or if it varies randomly with time due to bunching, as 
can easily be understood. The variation of shot noise 
intensity due to a randomly-varying event rate is similar 
to that of steady shot noise subjected to a random am¬ 
plitude modulation. The effect of amplitude modulation 
(AM) of a sinusoidal carrier is very well known to ra¬ 
dio engineers: modulation sidebands are produced below 
and above the carrier frequency. Mathematically, a car¬ 
rier at frequency v c that is AM-modulated at frequency 
v m develops sidebands at u± = v c ± z/ m : 

cos(27 iv c t) [1 + acos(27rz/ m t)] = cos(2-irv c t) 

+ ^ C0s(27TZZ + f) + ^ COS(27 TV_t) . (1) 

This result is easily generalized to a Fourier superposi¬ 
tion of modulation frequencies. Thus, the bunching noise 
component illustrated in Figure [3] may be interpreted as 
the sidebands on a DC carrier - the mean photocurrent 
- that are produced by the random modulation of the 
photon arrival rate. Indeed, these sidebands extend out 
to ±Azz which corresponds to the bandwidth of the ar¬ 
rival rate (or light intensity) fluctuations. Similarly, a 
Fourier component of shot noise at some frequency zq 
will develop sidebands extending over zq ± Azz as a re¬ 
sult of the arrival rate fluctuations. Because all Fourier 
components of shot noise develop these sidebands in the 
same way, it is clear that the net result must be a white 
spectrum. However, the sideband generation process in¬ 
troduces the possibility that the Fourier components of 
shot noise at different frequencies are correlated. Shot 
noise with a variable event rate arises in other contexts 
and is well studied, e.g . in the theory of diode mix¬ 


ers (Held & Kerr 1978) or in the dete ction of fast op¬ 
tical pulse trains (Quinlan et al. 2013). For these ex¬ 
amples, the event rate modulation is deterministic and 
periodic and the correlations between different Fourier 
components of shot noise play an essential role. However, 
these correlations vanish for the present case of photon 
bunching because the event rate varies randomly rather 
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than deterministically. Indeed, the photon shot noise 
must be a stationary process, in the sense that all statis¬ 
tics such as the autocorrelation function (J(f)J(f + r)) 
are invariant under time translation t —> t + t\, be¬ 
cause the time-varying photon arrival rate is also a sta¬ 
tionary random process. A translation in time changes 
the phase of a product of two different Fourier compo¬ 
nents, I*(y)i(v') -A r(y)i{v')e- i2 ^ v - v ' s > t \ and there¬ 
fore time-translation symmetry requires that the corre¬ 
lation of different Fourier components vanish. A more de¬ 
tailed mathematical proof of these sta tem ents is given in 

, which agrees 


section |dj culminating with e quatio n (j55][ 
w ith previous work (Picinbono, BenjdaBa 


_^ , , _ allah & Pouget 

1970). Thus, the fact that the shot noise spectrum is 


white tells us nothing about possible time-dependent 
variations of the shot noise intensity. So where is the 
bunching noise hiding ? To find it, we must go beyond 
the noise spectrum, which relates to the second order 
statistics of the photocurrent, and look at the fourth- 
order photocurrent statistics that are needed to describe 
the fluctuations of the shot noise intensity. This calcu¬ 
lation is presented in section [5] leading to equation (75), 
which reveals the presence ofthe bunching noise in the 
shot noise intensity. 


2. PHOTOCURRENT SPECTRUM FOR A SINGLE 
DETECTOR 

I now turn to a straightforward quantum-mechanical 
calculation of the output noise spectrum of a detector 
that is illuminated with filtered thermal radiation with 
bandwidth Av and occupation number n, as in the setup 
shown in Figure [2j The treatm ent uses a conven tional 
quantum formalism described in Zmuidzinas (2003); the 
calculations presented in this section are fairly standard 
and mainly serve to introduce the formalism and nota¬ 
tion. The principal result, stated below in equation (20) 
and illustrated in Figure |3j shows that the spectrum con¬ 
sists of three components: (1) a DC term corresponding 
to the average output; (2) a component due to photon 
bunching that is confined to a bandwidth equal to the 
optical bandwidth Au; and (3) a white-noise component 
due to photon shot noise that is limited only by the de¬ 
tector output bandwidth v,i- For conventional photon 
counting, the observable quantity is the time integral of 
the photocurrent Q which makes use of the fact that the 
DC photocurrent is proportional to the mean photon rate 
F = nAv. For the alternative shot noise measurement 
technique, the observable is the time integral of the noise 
intensity in the white-noise region Av < \v\ < z^, and 
makes use of the proportionality of the shot noise inten¬ 
sity to the mean photon rate T. 

Consider an ideal photon detector illuminated by a sin¬ 
gle mode of the radiation field. The radiation field is 


4 The time integral of the photocurrent is a useful and analyt¬ 
ically tractable quantity for quantifying the performance of ideal 
detectors but is not the optimal statistic for real detectors that 
have additional non-ideal sources of noise, e.g. amplifier noise. 
For example, one might use Wiener filtering followed by peak de¬ 
tection to locate and count the individual photon pulses in the 
output timestream; this approach rejects most of the noise ema¬ 
nating from the detection system during the time intervals between 
photon events. 


described by photon creation and destruction operators, 

poo 

&t(i)= / du e ~ i2nvt (jy) 

Jo 

poo 

b(t)= / dv e +i27VUt b{u) 

Jo 


which are defined only for positive frequencies and obey 
Bosonic commutation relations 

[b{v),tf{v f )\=6{v-v f ) . (2) 

I assume that the radiation field is in a thermal state 
described by the density matrix 


p = exp (^J dv j— x{v) b\v)b(v) + In 1 —e 


. ( 3 ) 

This is of the standard form for thermal equilibrium, 
p oc exp (—H/kT), where the Hamiltonian consists of a 
sum of harmonic oscillators, 


H = 


dvhv b\v)b(v) , 


( 4 ) 


and where x(v) = hv/kT is the normalized inverse tem¬ 
perature. The ln[l — e ~ x (")] term provides th e required 


normalization Trp =1. It is readily shown (Zmuidzinas 


2003) that this density matrix gives expectation values 
of 


(bJ {v)b(v')) = Tr \ptf {v)b(y')\ = n{v)8[y — v') (5) 

where the occupation number is given by the Bose- 
Einstein formula 


n{v) 


1 

e x W — 1 


1 

ghvjkT ^ 


( 6 ) 


Note that the excitation temperature need not be the 
same for all frequencies; we may easily generalize to 
x[y) = his/kT(v). 

An ideal photodetector produces one pulse at its 
output for every photon absorbed. For example, in 
a superconductor- insulator-superconduct or (SIS) tunnel 
junction detector (Tucker & Millea 19781, each absorbed 
photon causes one electron to tunnel across the junction. 
Such a detector may be described by a Hermitian pho¬ 
tocurrent operator 


Ip(t) = f dt'F{t — t')b\t')b{t') 

J — OO 


( 7 ) 


where F(t) describes the shape of the current pulse pro¬ 
duced by one photon. The detector output need not be 
an electrical current. More generally, we can consider 
/f(£) to be the output signal of the detector when illu¬ 
minated by the radiation field, and F(t) to be the out¬ 
put signal produced when a single photon is absorbed; 
however I will continue to call /p(f) the photocurrent 
operator. Note that I am assuming that the detector is 
operating in a linear regime: doubling the photon ab¬ 
sorption rate doubles the output signal. If the detector 
response is fast relative to the timescales of interest, we 
may approximate 

m « s(t) , ( 8 ) 
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and therefore consider the Hermitian operator 

m = b'mt) ( 9 ) 

where now I(t ) has units of s _1 . 

The impact of photon bunching on the sensitivity 
of measurements performed using conventional photon 
counting may be demonstrated by considering the oper¬ 
ator 



which represents the number of photons detected in a 
measurement time interval [— T/2,T/2], The mean and 
varia nce of this operato r are given in equations 41 and 
42 of Zmuidzinas (2003): 


(ID 


POO 

(Nt) = Tr (pNt) = T / dvn(v) 

Jo 

pOO 

^Nt = ( n t) - (Nt ) 2 = T / dvn(v) (1 + n(v)) . 

Jo 


Thus, the fractional measurement uncertainty is 


cr Nt \/l T fi 

m ~ 7m 


(13) 


and is degraded by y/1 + n due to photon bunching 
as compared to the fractional Poisson uncertainty of 
l/\/ {Nt)- Here the effective occupation number h is 
defined by 


- fo°° dl/ [niy)f 

Jo°° dvn{v) 


(14) 


For the simple case that n(v) = n inside an optical band¬ 
width A v and zero outside, one readily finds n = n. 

The calculation of the output noise spectrum of the 
detector makes use of the Fourier transform of the pho¬ 
tocurrent operator: 


i{v) = P{-v) 

/*+ oo 


dt e~ l2nvt I{t) 


' —OO 

/*+ oo 


dte 


poo 

dv i / dv2b\vi)b(v2) 
Jo 


to Jo 

/*+oo 


dj. e +i2-KVit ^—i2'KV2t ^—i2'Kvt 


p oo 

/ dv i tf(v i + v)b(vi)0(vi + v) . 
Jo 


(15) 


The photon operators b and are only defined for posi¬ 
tive frequencies; the unit step function 9(vi+v) is needed 
to guarantee v\ + v > 0 even when v < 0. For brevity of 
notation, I will not write the step function explicitly but 
instead rely on the interpretation b{y\) —> b(vi)6(vi) and 
similarly for Clearly, the effect of a finite pulse 


width F(t) will be multiplication by the corresponding 
frequency-domain filter F(y ), 


I F (v) 



dte- anvt I F (t) =F{v)I{v) . 


(16) 


The power spectrum Si(v) = Si(—v) of the photocurrent 
is dehned by 

(P(v)i(v')) = S^v) d(v - v') ; (17) 

including the pulse shape F(t) would lead to a filtered 
power spectrum 

S If (v) = \F(v)\ 2 S i (v) . (18) 

Because the thermal density matrix is gaussian, the re¬ 
quired expectation value may be found by combining the 
photon operators pairwise, 

dvidv2 

(5 t (^i)&(^i + v)b\v 2 + v')b(v 2 )) 
dvidv 2 

(b\vi)b(yi + v)) (b ] (v 2 + v')b{v 2 )) 

+ (tf(vi)b(v 2 )) (b(v i + v)tf{v 2 + v')) 

POO 

= S(v — v') / dvidv 2 {n(vi)n(v 2 )S(v) 

Jo 

+ n(v l )[n(vi +v)+ \}5{vi~ v 2 )} . (19) 

The photocurrent power spectrum is therefore given by 
a sum of three terms, 





Si{v) = T 2 5(v) + r + / dvi n(vi)n(vi + v) , (20) 


where, according to eqn. (11), the mean photon arrival 
rate is 


r = 


dv i n(vi) . 


( 21 ) 


When n(v) is constant within a bandpass Av and zero 
outside, we obtain T = nAv, so the occupation number n 
may be interpreted as the number of photons per second 
per Hertz of optical spectrum. 

The first term in Sj(v ), proportional to 5{y), repre¬ 
sents the contribution to the power spectrum from the 
DC value of the photocurrent. At nonzero frequencies, 
only the second and third term contribute. The second 
term, F, is white noise independent of frequency v , and 
represents photon shot noise. The third term is due to 
photon bunching, and is not white. Indeed, for a rectan¬ 
gular optical bandpass of width Av, the spectrum of the 
bunching term has a triangular shape that is symmetric 
with respect to v = 0 and extends over — Av < v < +Av: 


^(bunching) 


\ nT 

i M 

1 Av 

l 

0 


\v\ < Av 
otherwise 


( 22 ) 


The sum of these three terms is plotted in Figure [3] At 
high frequencies \v\ > Av, only the white photon shot 
noise term contributes. This suggests the following idea: 
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place a high-pass filter F{v) at the detector output that 
transmits only at frequencies \v\ > Av. The spectral 
density of the shot noise is F; we can therefore measure 
the photon rate by measuring the noise intensity at \v\ > 
Av. However, as we shall see, this method does not avoid 
the sensitivity degradation due to photon bunching. 

3. PHOTOCURRENT CROSS-SPECTRUM FOR MULTIPLE 
DETECTORS 


thermal 

source 

Av 



termination 


T « hv/k 90° hybrid 
coupler 


~>h (t) 


detector 1 


-((d)— >h(t) 

detector 2 


Fig. 6.— The radio-frequency equivalent of the two-detector plus 
beamsplitter setup shown in Figure[l] Free-space propagation is re¬ 
placed by guided-wave propagation in transmission lines or waveg¬ 
uides, the funct ion of the be amsplitter is performed by a 90° 3 dB 
hybrid coupler < |Pozar 2012), and the unused port of the c oup ler is 
connected to a cold termination. As shown by equation ( |41[ ), the 
spectrum of the difference of the photocurrents I/± ( t ) = I\ (t)—l 2 ( t ) 
does not contain the DC or bunching noise components, and is 
therefore white across the full detector output bandwidth. 


I now generalize the discussion of section [2] to the case 
of multiple detectors, in order to analyze the two-detector 
scheme illustrated in Figure [l or its radio-frequency 
equivale nt shown in Figure [b| T le detection scheme pro¬ 
posed by|Lieu et al. (2015 [) uses an identical two-detector 
setup. I start with a more general case in which an arbi¬ 
trary passive linear optical system is used to illuminate a 
set of detectors, and calculate the cross-spectral density 
of the output currents. The principal result (equa tion 311 
is a straightforward generalization of equation (201 tor 


a single detecto r; to ourknowledge this is a new re¬ 
sult. Following Zmuidzinasl (2003), the optical system 


is represented by a passive linear 7V-port network with 
a scattering matrix S; a 50/50 beamsplitter is an exam¬ 
ple of a four-port network. The network is illuminated 
by incoming radiation described by the photon operators 
ai{y), and produces outgoing radiation according to the 
scattering equation 


bi{v) = ^2 Sij(v)ai{y) + d(v) 

3 


(23) 


where the d{v) are operators representing noise added 
by the network. Here the indices 1 < i,j < N label 
the ports of the network. The noise operators satisfy 
commutation relations 

{d(v),Cj{v')\ = [1 - S(v)S\p)]jiS{v - v') , (24) 

as required to preserve the Bosonic commutation rela¬ 
tions for the output operators, 

[bi(v), b]{v')\ = 6 i:j 5(u - v') , (25) 

given that the input operators also satisfy the same com¬ 
mutation relations, 

[ai(v),a](v')]= 5ij5(v - v') . (26) 


If the input radiation is thermal and the TV-port is 
passive, the output radiation is also thermal and may 
be fully described by a mode occupation matrix Bij(y ), 
defined through 

(&i(!z)&j(z/)) = Bij(y)5(v - v') , (27) 

which is a generalization of the mode occupation number 
n{v). Photocurrent operators and their Fourier trans¬ 
forms may be introduced for each port: 

I i {t) = b\(t)b i (t) . (28) 

nOO 

ii(y)= / dvi b\{v-i +v)bi(i/i) . (29) 

J o 

The photocurrent power cross-spectrum C t j (y) is defined 
through the expression 

= Cij(v)5{v - v') , (30) 


which I calculate using pairwise evaluation of the result¬ 
ing fo urt h-order moments of the photon operators, as for 
eqn. (191. The power cross-spectrum is found to be 


Cij(v) = r^- 5{v) + Tj Sij 

/»oo 

+ / dui Bij(yi)Bji(y\ + v) , (31) 

Jo 


where the mean photon rates at the detectors are 


r, = 


nOO 

/ dv 1 B ii (v 1 ) . 
Jo 


(32) 


As was found for the single detector case and illustrated 
in Figure pi we see that the power cross-spectrum Cij(v) 
for multiple detectors consists of three terms: a DC term 
(f itjSiy)), a white spectrum from photon shot noise 
that is uncorrelated between detectors (f<5,; 7 ), and a 
bunching spectrum 

/»00 

^(bunching) {v)= / dviB .. {yi)B .. {yi + v) (33 ) 

Jo 

that exhibits correlations between detectors but is lim¬ 
ited to frequencies v < Av. 

We now specialize to a four-port netwo rk a ppropriate 
for a beamsplitter or 90° 3 dB coupler (Pozar 2012), as 
illustrated in Figures [l] and [6j with a scattering matrix 
given by 

0 0 * 1 ' 

OOli 
i 1 0 0 
1*00 


\/2 


(34) 


This matrix is reciprocal, S T = S, as required by time- 
reversal symmetry, and S is also unitary, SS' = 1, and 
therefore the network is lossless. We will call ports 1 
and 2 the output ports and place detectors on them, and 
ports 3 and 4 will serve as the input ports. The incoming 
fields aiiy) are assumed to be in independent thermal 
states described by occupation numbers rii(u). Port 4 
will be illuminated with occupation number ni(i'). Port 
3 will be terminated in a vacuum state with zero occupa¬ 
tion number; furthermore, the detectors are assumed to 




























be cold and therefore do not radiate toward the beam¬ 
splitter, so n\{v) = ri 2 {v) = 773 ( 7 /) = 0. We may now 
calculate 


= ^ s ik{v)S* k {y)n k {v) . 
k 



= S i 4 ,(y)S* i (y)n i (v) . 

(35) 

Therefore, 


Bn {v) = B 22 {v) = UA ft 

(36) 

while 


Bn{u) = B* 21 {v) = -i^ . 

(37) 


Thus, the bunching power cross-spectrum for the detec¬ 
tor ports 1, 2 is 


^(bunching) _ _ 

~ 4 


i i 
i i 


disi 714 (^ 1 ) 714 ( 1/1 + v) ; (38) 


which reproduces an ear lier result of |Picinbono, Benjda-| 
ballah & Pouget| (fl970|; the purpose of presenting the 
detailed derivation here is to introduce the formalism 
in preparation for the calculation of shot noise intensity 
fluctuations in section 0 

A classical current containing shot noise, e.g. the cur¬ 
rent across a tunnel barrier with a low transmission prob¬ 
ability, may be considered to be a sum of impulses, 


!{t) 


(42) 


where {tj} represent the times at which discrete charges 
(e.g., electrons) jump across the barrier. To remain con¬ 
sistent with sections [2] and [3j I have omitted the usual 
factor of electron charge e, so I(t) has units of s -1 or 
Hz. Suppose further that the average current is time- 
dependent, 

(I(t))=T(t)=f + ST(t) , (43) 


the bunching noise is fully correlated between the detec¬ 
tors. The matrix in this expression has eigenvalues of 2 
and 0 for the symmetric and antisymetric eigenvectors 
(1,1) and (1,-1). We therefore see that the bunching 
noise term will be absent for the difference of the two 
detector photocurrents, I a (t) = h(t) —12 (t). Neglecting 
the DC term, the noise matrix is 


C(*) = 


1 

2 


1 0 
0 1 


du 1774 ( 7 / 4 ) 


+ 


1 

4 


1 1 
1 1 


dv 177 . 4 (^ 1 ) 714(^1 + v) 


(39) 


where the first term represents the shot noise, which is 
uncorrelated between detectors. This is the only term 
that remains when we calculate the spectral density of 
/a, 




[ 1 -1 ] C(u) [ 

nOO 

/ dv\ 774 ( 7 / 4 ) = f , 


(40) 


where P is the mean event rate and <5T(f) represents vari¬ 
ations in the rate and therefore has zero mean. If the 
event rate T(f) is constant, i.e. ST(t) = 0, we know 
that the current has a shot noise spectrum that is white 
and has intensity T. On the other hand, if the event 
rate varies with time, we expect the shot noise intensity 
to also vary. Thus, a time-resolved measurement of the 
shot noise intensity should allow us to measure the cor¬ 
responding time-dependent current. This possibility is 
investigated here and furth er in section [5| 

I follow the approach of Kelley & Kleiner (1964) for 


calculating classical shot-noise statistics; a simil ar but 
mathem atically m ore for mal a pproa ch is given by[Picin- 


bono, Benjdaballah & Pouget (1970). Let y, be a random 
variable that represents the number of charges that flow 
during the time interval [ti,U + At,]. I assume that yi 
is independent of all other yj for j 7 ^ i, and furthermore 
that for sufficiently small time intervals At,. yi has a 
probability distribution given by 


P(yi = 1 ) = T(ti)Ati 
P(yi = 0 ) = 1 - T(ti)Ati . 


which is the same as the shot noise intensity for a single 
detector without the beamsplitter. To summarize, the 
difference of the two detector currents / a (£) = h (t) — 
/ 2 (f) has a white spectrum, given by 

(li{v)iA{v'))=T 5 {v- v ') , (41) 


The number of charges that cross during the time interval 
[0, T] is a random variable given by 

pT M 

N t = I{t)dt^Y^ Vl (44) 

i=1 


and illustrated in Figure [4j 


4. SPECTRUM OF SHOT NOISE WITH A VARIABLE 
EVENT RATE 

As discussed in section [lj in the semiclassical picture 
one views photon bunching as being caused by a stochas¬ 
tically varying photon arrival rate. Here I calculate the 
spectrum of classical shot noise for a time-varying event 
rate and demonstrate that the shot noise remains white 
and uncorrelated, as discussed qualitatively in section [l 
I further demonstrate that the output noise spectrum of 
a single detector calculated in the semiclassical picture 
can reproduce the quantum-mechanical result give n i n 
section [2j The principal result is given by equation (55), 


where the time interval [0, T] has been split into M 
nonoverlapping subintervals [£,;, t, + Aij]. The distribu¬ 
tion of Nt is encoded by the moment generating function 
calculated in Appendix |A} 


M 


G N (s) = (e sNT ) = J2h 51 {Vix-Vin) 


k\ 

k =0 i 


= exp \y T (e s - 1)] 

t4 

k! 


(45) 


00 h 

ST fA. 
j,i c 


k =0 


As expected (Mandel 1958), this result shows that Nt 
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follows a Poisson distribution with mean 

r T 


Ht = f dtT(t) . 
Jo 


(46) 


Thus, the current /(f) is a Poisson process with a time- 
dependent rate F(t). 

The same formalism can be used to work out the spec¬ 
trum of shot noise for a time-dependent current. The 
time-limited Fourier transform of the current is defined 
as 

r + T /2 


It{v) = 


l-T/2 


dtl{t)e 


—i2ivvt 


(47) 


which allows the power spectrum to be computed by eval¬ 
uating the limit 


lim (It(v)I^(v') 
t— x- ' 


(48) 


Expressing the current in terms of the random variables 
Vi gives 


M 


i T {v )« ^2 Vi?' 

i—1 


Therefore, 

[i T {v)W) 

Now 


M 


Yum) 


—i2-Kvti -\-i2 , Ki''tj 


(49) 


(50) 


*.j=l 


{viVj) = (Vi) (vj) + s ij (ivi) ~ (yiY 


+ 5i 


r(t ! )At I -r 2 (t i )(At l )" 


(51) 


making use of y 2 = yt. Inserting the first term into the 
sum and taking the continuum limit gives 


lim 

T—> oo 


,T/2 r T/2 

/ dh r(ii) e - i27ri/tl / dt 2 T{t 1 )e +i2 ^ ti 
J-T/2 J-T/2 


= (r<y(i/) + rf^i/)) (rs(i/) + 8f * (i/)) (52) 

while the second term yields 

T8(v — u') + ST(v — v') ; (53) 

note that the (A tj) 2 term vanishes in the continuum 
limit. Thus, the Fourier components of shot noise are 
correlated when the event rate varies deterministically 
with time, according to 

i(v)i* (i/)) = (r$(i/) + st (o) (f«(i/) + tff*(i/)) 


+ T8{v — v') + <5r(^ — v') . 


(54) 


If the event rate has stochastic time-dependent rate 
variations <5F(t), taken to be a stationary random process 
with power spectrum Sriy), the sho t no ise spectrum may 
be obtained by averaging equation (54) over JT: 

i{v)t{v')\ = [f 2 ^)+f + S r (v)\ 

/ y,sr 

x 8{v — v') ; (55) 


the 8{y—v') factor indicates that different Fourier compo¬ 
nents are uncorrelated and therefore this compound Pois¬ 
son process is stationary, as expected. Specifically, the 
shot noise component T5(;/ — v') remains whi te a nd un¬ 
correlated, as promised in secti o n[l| Equati on (55) agrees 
with the result of Picinbono, Benjdaballah "SPPouget] 
(1970) (their equation 2.27), who claim agreement with 
an earlier result by Mandeb 
The power spectrum of this classical compound Pois¬ 
son process may be made identical to the spectrum of the 
pho tocu rrent calculated quantum-mechanically (equa¬ 
tion 20), provided we make the identifications 


p+OO 


r = 


dv\n(yi) 


(56) 


and 


r+oo 

S'r(^) = / dv\n(v\)n{y\ + v) . 

Jo 


(57) 


If the occupation number n{v) is constant across an op¬ 
tical bandwidth Av and zero outside, the relative im¬ 
portance of the bunching and Poisson terms at noise fre¬ 
quencies well below Av is governed by 


Sr(0) _ C°° dvm 2 ^) 


/ 0 +o ° dvin{v\) nAv 


(58) 


5. SHOT NOISE MEASUREMENT: SEMICLASSICAL 
ANALYSIS FOR A SINGLE DETECTOR 

I now turn to the computation of the fluctuations in 
the intensity of classical shot noise with a time-varying 
event rate and demonstrate that the shot noise intensity 
reflects variations in the event rate. I apply these results 
to the case of photon detection under the assumption 
that bunching may be described by a stochastic pho¬ 
ton arrival rate whose mea n and po wer spectrum are 
described by equations ( [56]) and ( [57]) , respectively. As 
discussed in section [l] and Appendix |F] this assumption 
is consistent with the full quantum theory and is sup¬ 
ported by experiment. T he principal resu lts presented 
in this section (equations [68 71 72 and 75) are new 


and demonstrate that shot noise measurements using a 
single detector suffer the same y/1 + n sensitivity degra¬ 
dation due to photon bunching as would occur for direct 
photon counting. 



detector 


Fig. 7.— Signal flow diagram for shot-noise detection using a sin¬ 
gle detector. A noise filter W(v) allows selection of the portion of 
the spectrum where shot noise dominates (e.g., the hatched region 
in Figure: fui prior to the measurement of the noise intensity using 
a square-law detector and integrator. 


Because the rate fluctuations Sr(v) have a limited 
bandwidth (equation [57|) , only shot noise contributes to 
the current noise (equation 55) at high frequencies, and 
therefore a measurement of the intensity of the high- 
frequency noise should give us information on the mean 
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rate P and therefore the mean current. Specifically, the 
output noise spectrum for a photon detector as illus¬ 
trated in Figure [3] suggests use of the shot noise mea¬ 
surement setup shown in Figure[TJ which includes a noise 
filter W(v) for isolating the shot noise dominated portion 
of the spectrum. We therefore consider applying a filter 
to the current, 


However, we expect that the intensity of the shot noise 
should be affected by the rate fluctuations 5T(t). We are 
therefore interested in the fluctuations 

a p = (Pw,t) ~ ( p w,T ) • (64) 

This quantity will require evaluation of the fourth mo¬ 
ments of the form 


Iw{t)= f , 

J —oo 


(59) 


and then continuously integrating the noise power over 
a measurement interval [—T/2, T/2], 


Pw,t = 


,T/2 

/—T/2 
,T/2 

/-T/2 

,T/2 

/-T/2 
/*+oo 


dt [Iw(t )f 


dt 


dt 


dt±W(t — t\)I(ti) 


\_J —OO 
/*+oo 


dve l27rvt W{v))I{v) 


dv i dv 2 I (vi)I* (v 2 )W (vx)W* (v 2 ) 


,T/2 
'-T/2 

/* + OO 


dt e i2nVlt e~ i2 ' KV2t 


(60) 


dm 


i(^i) |VF(i/i)| as T — > oo. 


The mean value of this measure of the high-frequency 
noise intensity is 


( P\V,T) = 


f-\-oo 


+oo 


dv i 


f T/2 

-T/2 

/»+oo 




dv^s^iw^Y 


dv 2 (i{v 1 )i*(v 2 ))W(v 1 )W*(v 2 ) 


r T/2 


dt 


I-T/2 


p+oo 


= rr 

+ T / 




^ —OO 

+oo 


di/i |W(i^i)| S'r(^i) . 


( 61 ) 


The second term picks up the bunching noise component 
but may be made negligible by choosing W(y) to be zero 
at the lower frequencies where Sr(v) has an appreciable 
value, as illustrated by the hatched region in Figure [3j 
With this choice, if the noise measurement bandwidth is 
defined as 

/• + 00 


2 B = 


dv ! \W(yi)\ , 


(62) 


where the factor of two accounts for negative frequencies, 
we have 

(P W , T ) = 2TBT , (63) 

and therefore the measurement scheme shownhr Figure [T] 
provides us with the desired information on T. 


(i{vr)i{v 2 yi{v3)i{v A y) . (65) 


Because of the presence of the high-pass filter W{v), we 
may safely assume that none of the frequencies are zero, 
and therefore omit the DC terms. Using the same ap¬ 
proach as before, we write 


(i{vi)i{v2)*I{v?,)i{vi)*y ~ {ViViVkyi) 

V ijkl 

x ^-12-Ki/^ti e +i2^f 2 tj e -i2-Kf 3 t k ^-12-KViti j-gg-j 

Appendix [5] provides the details of the evaluation of 
this quantity, leading to an expression for <r 2 P in the 
long measurement time limit A vT >> 1 involving seven 
terms, labeled Alb+c, B2+B3+B4-I-B5, C2a+C3a, Clb, 
C2b+C3b, C4b+C5b-|-C6b-|-C7b, and Dl. Three of 
these terms drop out if we design our filter W ( v ) so that 
it rejects noise due to the rate fluctuations, i.e. 



dv |W(^)| 2 Sr{v) —> 0 , 


(67) 


as illustrated by the hatched region in Figure[3| The sur¬ 
viving terms (C2a+C3a, Clb, C2b+C3b, and Dl) con¬ 
tribute fractional fluctuations of 


4 _ 1 f 2fdv\W(v)\ 4 Sr_(0) 

( Pw,T ) 2 ^ \ [/ dv\w(v)\ 2 ] 2 r 2 

2 / dvdv'\W(v)\ 2 \W(v')\ 2 S r (v - v') 1 

f 2 [f dv\w(v)\ 2 ] 2 r 


( 68 ) 


The interpretation of these terms is simplified by choos¬ 
ing a filter function W ( v ) which is unity inside a mea¬ 
surement bandwidth B and zero outside, so that 



dv\W(v) | 2 



dv\W(v)\ 4 = 2 B , 


(69) 


including contributions from positive and negative fre¬ 
quencies. If we define an effective bandwidth for the rate 
fluctuations, 


f dvS r (v) 

S r (0) 


(70) 


and evaluate the third term under the assumption that a 
wide bandwidth is chosen in order to optimize the shot- 
noise measurement, B >> Av , the terms simplify to 


a% = JL , M 5 ! + g r(0)A^ 

(P WT ) 2 BT + T 2 T T 2 BT + TT ‘ 


(71) 


The last term is due to the Poisson fluctuations in the 
number of events over a time T that one must have even 
if the event rate is constant. Meanwhile, the first term 
represents the noise that results from measurement of 
a finite number of independent samples associated with 
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the time-bandwidth product BT\ indeed, one sees that 


this term reproduces Dicke’s result (Dicke 1946|) in the 
shot noise context. Thus, it is helpful to use a large 
shot noise measurement bandwidth, although the Pois¬ 
son term dominates when the bandwidth exceeds the 
mean event rate, B > T. The second and third term 
represent the effect of event rate fluctuations, being pro¬ 
portional to the spectral density of the fractional fluctu¬ 
ations Sr(0)/T 2 . Here the spectral density Sr(v) is eval¬ 
uated at zero frequency because the rate fluctuations are 
being averaged over a long meas urement time T. The 
last three terms in equation © may be rearranged to 
read 


1 


(Pw,i 


1 + 


S r (0) S r (0) Az/ 

f + f ~B 


■ (72) 


In this form, the last term in the square brackets can be 
seen to represent a correction to the rate fluctuation term 
due to finite measurement bandwidth and is negligible 
under the assumption B » Av. 

If we make use of the identifications a ppro priate for 
thermal photon noise given by equations (56) and (57), 
and furthermore assume that the occupation number 
n{v) is constant inside an optical bandwidth Av and zero 
outside, we have 

f = nAv (73) 


and 


S r (0) = n 2 Av . 


(74) 


The fractional fluctuation in the shot noise intensity is 
then given by 


(Tp 


( Pw,T) 


1 

BT 


1 + n 

Tt 


(75) 


In t he l imit B » T, we recover the usual result (equa¬ 
tion 13) that photon bunching gives a sensitivity penalty 


of y/1 + n as compared to Poisson statistics. This occurs 
despite the use of the white portion of the shot noise 
spectrum to measure the photon rate. 


6 . 


SEMICLASSICAL ANALYSIS FOR MULTIPLE 
DETECTORS 


The extension of the semiclassical treatment in sec¬ 
tion [5] to the case of multiple detectors is straightforward 
and allows us to analyze the sensitivity of shot noise mea¬ 
surement schemes ap pli ed t o the two-detector setup pro¬ 
posed by |Lieu et al.| (|2015|), shown in Figures fl] and [6| 
The principal results (equations 90 and 93) are new ana 
agree with those in section [5j theyaemonstrate that shot 
noise measurements appliea to the two-detector scheme 
also cannot evade the y/1 + n sensitivity degradation due 
to photon bunching. We perform our analysis for a signal 
processing setup (Fig. [8]) similar to those ty pically used 
f or experim ental measurements of shot noise (Schoelkopf 
et al. 19971, though i t differs in detail f rom the signal pro- 


cessmg proposed by Lieu et al. (2015). Nonetheless, our 
calculations are directly applicable to the regime that 
Lieu et al. claim leads to suppression of the bunching 
noise; a detailed comparison of the calculations is given 
in section [HI 


Suppose we have multiple currents exhibiting shot 
noise, 

Ig(t) = y ^6{t- ti >a ) (76) 


with time-dependent event rates 

(i a (t)) = r a (t) = r a + sr a (t). 


(77) 


Here the rate fluctuations are stationary stochastic pro¬ 
cesses described by a cross-spectral correlation matrix, 

^(^(z/)) = C£\u)8(v - v') . (78) 

Here a and b are discrete indices that label the currents. 
As before, I discretize time and introduce random vari¬ 
ables y ai i to represent the number of events for current a 
in the time interval [t*, At*]. The cross-spectral density 
between two currents is given by 

Ial(vi)ia2(v2)) ~ 




x e -27r!yiti e +27ri/ 2 t i 


The y a ,i are all independent, so 

(iJaA XUaM.j) = ( Val,i ) {Va2j) 

+ ((j/al.i) — ( Val,i ) 


(79) 


(80) 


The term ( y a i,i ) 2 is of higher order in AL and can be 
neglected in the continuum limit: 

= 


(81) 


(r al (i(zq) + ( 5 f Q i(zq)) (r a2 <5(z/ 2 ) + ^*2(^2)) 


+ Sal, 0,2 r a ici(zq - V 2 ) + 5fal(^l “ ^ 2 ) 

Averaging over the random process i5r a (t) yields 
7al(l / i)/* 2 (z'2)') = [f alfa 2 ^(^l) + <ial,a 2 ? al 

/ y,SF 

+ C 'al,o2( I/ l) - V 2 ) , (82) 

which is simply a gener aliza tion of the single-detector 
result given in equation (55). If we compare this result 
to equation (31) for the photocurrent correlations among 
detectors illuminated with thermal radiation, we see that 
the expressions coincide if we make the identifications 


T a = I a = / dvB aa (y) (83) 

Jo 


and 


/•OO 

C ab\ v ) = / dv B ab (y')B ba (y' + v) 
Jo 


(84) 


which are generalizations of equations (56) and ( |57| ) . 

The shot-noise measurement scheme fora singleaetec- 
tor shown in Figure[7]may easily be adapted for use with 
two detectors as shown in Figure |8j this setup is designed 
to measure the shot noise intensity in the difference of 
the two currents, I a = h — I 2 , as proposed by [Li eu et] 


al. (2015). Although the filter W(z') is no longer needed 
Tor rejection of the bunching noise at low frequencies, it 
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A(t) 





--Pa 


filters differencer 

Fig. 8.— Signal flow diagram for shot-noise detection using two 
detectors. In principle, the ordering of the noise filtering and dif¬ 
ferencing operations may be interchanged since both are linear. 


is maintained in the setup since any real system has a 
finite bandwidth. The output of the shot noise intensity 
measurement is given by 


I dti ( dtiW(t-ti)[Ji(ii)-J 2 (ii)]} 

J—T/2 U-oo J 


r T/2 

Pa = I dt- 

l-T/2 

f+OO 

dis 1 dis 2 W ( v\) W* (v 2 ) 


h{vi)It{vi) + h{vi)I*2 (4 - 2/ 1 (iv 1 )I 2 > 2 ) 

rT/2 

x / dte i2 ™ lt e~ i2 ™ 2t . 


J-T/2 

and has an average value 

/*+oo 


(85) 


(Pa)=T I dv\W(v)\ 2 T 1 + C$(is) + T 2 


+Cg\v)-2C'&(v) 


( r )f 


( 86 ) 


Equations ( |38[ ) and ( |84[ ) give 

') = 

dis' 714 ( 1 /) 77 . 4 ( 1 / + is) 


c£\v) = Cg\u) = C&\u) 


\ 


and therefore 


(P A ) = 2 BT (f! + f 2 ) 


(87) 


( 88 ) 


is a measure of the total event rate regardless of the 
choice of the filter passband W (is). 

Calculation of the sensitivity of this shot-noise inten¬ 
sity measurement requires evaluation of fourth-order mo¬ 
ments of the photocurrent, 

F abed = (l™I< w >lWl< w >) 

= f dis 1 dv 2 dv3dis±h a (is 1 )IZ(v 2 )L(i'3)i*d(i'4) 

• 7—00 ' 

x W(is 1 )W*(p 2 )W(is 3 )W*(is a ) 

rT/2 

x / 

J-T/2 

r T / 2 

x / dt' e i2 ™ 3t e- i2 * Vit . (89) 

J-T/2 


This expression may be evaluated using the same ap¬ 
proach as used for the second moment; the details are 


given in Appendix [C] The resulting variance of the shot 
noise intensity is derived at the end of the Appendix: 

= (Pi) - 4a) 2 

= 2-Fiin 4- 2Fn 22 + 4 Fi 2:L2 — 8Fm 2 


= +8TT 2 / dis\W(is)\ 4 +4TC[ r 1 ) (0) 


dis\W(is)f 


+ 8 T J disdu'\W(is)\ 2 \W(is')\ 2 c[l\is - u') 

dis\W(is)\ 2 


+ 2TT 1 


(90) 


To compare to the previous single-detector case (equa¬ 
tion HD , we make the substitutions 


and 


C 


(r) 


Ti = -T 
2 


M = ^r(i/) ■ 


(91) 


(92) 


Using the mean value of Pa given by equation (86), we 


may express the fractional fluctuation in the noise inten¬ 
sity of the difference current as 


u Pa 

4a) 2 


i 

T 


2 f dis\W(is )| 4 S r (0) 


[Jdis\W(v)\< 


r 2 


2 f disdis'\W(is)\ 2 \W(is')\ 2 Sr(is - is') 1 

f 2 [J dis\W(is)\ 2 ] 2 r 


(93) 


which is exactly ou r pr evious result for a single detector 
given by equation (68). In particular, the rate fluctua¬ 


tion term S'r(O) leads to the \/l + n sensitivity degrada¬ 
tion due to photon bunching. However, it is no longer 
necessary to make the assumption that the filter W(is) 
rejects the low-frequency excess noise; differencing the 
two detectors fed by a 50/50 beamsplitter takes care of 
the rejection instead. 

7. SHOT NOISE MEASUREMENT: A QUANTUM 
CALCULATION 

The semiclassical analyses given in sections [5] and [6] 
for shot noise measurements are revisited in this section, 
but now making use of a fully quantum-mechanical treat¬ 
ment. I focus first on the single-detector case; the gener¬ 
alization to multiple detectors is straightforward and is 
given at the end of this section. As before, the photocur¬ 
rent operator is given by equation (19]) and has Fourier 
components given by equation <E1- ^ he shot noise in¬ 
tensity is measured in the same way: a filter W(is) is 
applied to the photocurrent before a square-law detector 
and integrator are used to measure the intensity, as illus¬ 
trated in Figure [7] This measurement sc hem e produces 
the quantity Pw,t as defined by equation (60). However, 
calculation of the statistics of Pw,t now requires quan¬ 
tum operator averages, which I perform in the usual way 
appropriate for thermal radiation, namely by combining 
photon creation and destruction operators pairwise. The 
quantum computat ion of the second-order moment is de¬ 
tailed in equation (19), with a result that is identical to 
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the semiclassical second-order moment, 

(7(^)1 t(i/ 2 )> = [f 2 5(n)+f + 5 r (^i)] 

x <y(^! - i/ 2 ) . (94) 

In fact the semiclassical spectrum was chosen to coincide 
with the quantum result t hrou gh th e de finitions of T and 
Sr(v) given in equations (561 and (571. Thus, we con¬ 


clude that ( Pw,t) = TBT for the quantum calculation 
just as for the semiclassical case (equation [60|, provided 
that we choose the noise filter W(is) to avola the excess 
low-frequency noise as discussed in section [5] 

Evaluation of the fluctuations of Pw,t requires a quan¬ 
tum computation of the fourth-order moment of Fourier 
components of the photocurrent, which in turn requires 
eighth-order moments of the photon operators: 

F(v\,V2, V3, v±) = (i/ 2 )/(i / 3)/ t (rq)) 

= J (lv\ dv'^dv'.^dv'^ 

x (v[)b(v[ + vi)b ] {v 2 + p 2 )&K) 

+ + . (95) 

Combining operators pairwise produces 4! = 24 terms. 
However, as for the semiclassical calculation, many of 
these represent DC terms that are rejected by the fil¬ 
ter W(v) and therefore do not contribute to the shot 
noise intensity. For example, if the first two operators 
are paired, we will have a factor 

i 11 ') = J du[ {b ] {v[)b{v[ + u{)) 


= J dv[n(v[)5(v{) 


= f5(i/i) 


(96) 


which vanishes except at DC, v\ = 0, and may there¬ 
fore be ignored. This shows that we may ignore any 
similar pairing, e.g. (22') in the obvious notation. 
Any pairing may be represented by a permutation, 
e.g. (11') (23') (34') (2'4) corresponds to the permutation 
(1)(234) expressed in cyclic notation. All permutations 
that include a cycle of length 1, e.g. (3), will give DC 
terms that we may ignore. This leaves 9 permutations 
left to consider: 


(12) (34) (1234) (1243) 

(13) (24) (1324) (1342) 

(14) (23) (1423) (1432) 


(97) 


The detailed evaluation of some of these pairings is 
given in Appendix[DJ These pairings reproduce the terms 
found in the semiclassical calculation outlined in sec¬ 
tion [5] and detailed in Appendix [B] and also generate 
some extra terms that arise from the non-gaussianity of 
the photon arrival rate fluctuations as described in Ap¬ 
pendix [F] that are neglected in our semiclassical calcula¬ 
tion. In particular, the (1432) permut ation inclu des the 
contributions expressed by equations (D5) and (D6) in 
Appendix |D| 


( P ^>(1432)= T [ f + ^(0)] [| dv\W(u)\ 


which are same as the semiclassical terms D1 and and 
Clb listed in Appendix [B] that correspond to Pois¬ 
son noise and bunching noise, respectively. The latter 
term contributes Sr (0)/T 2 to the fractional fluctuations 
Opf {Pwt) ( see equations 


68 


75) and thus represents 


the shot noise intensity fluctuations due to photon bunch¬ 
ing. It is this term that gives the same \/l + n sensitivity 
degradation due to bunching as occurs for ordinary pho¬ 
ton counting. 

It is not difficult to translate these results to the case 
of multiple detectors. We again focus our attention on 
the (1432) operator pairing in the corresponding quan¬ 
tum calcu lati on, w hich includes the contributions (equa¬ 
tions 


D7 and D 81 


[Fabcd] (1^32) — F^ab^cd 

+ ... 


'f„*ac + C£>( 0 ) 


1 2 


dv\W(v)\ 

(99) 


that correspond to the Poisson and bunching terms D1 
and Clb found in the semiclassical calculation, as out¬ 
lined in section[6]and detailed in Appendix |C| The latter 
term contributes 


ATCii (0) 


dv\W{y) | 2 


( 100 ) 


to the measurement variance Up A (equation|90j) for the 
two-detector setup shown in Figures [l] and [d|and leads 
to the vTTn photon bunching degradation. 

Thus, we conclude that a full quantum calculation re¬ 
produces the conclusions of the semiclassical analyses for 
one or two detectors given in sections [5] and [6j namely 
that shot noise intensity measurements cannot evade the 
bunching noise. 

8. CORRELATION OF SHOT NOISE AND PHOTON 
COUNTS 

We have two ways of measuring the photon flux: direct 
photon counting using the time inte gral of the photocur¬ 
rent, Nt, defined in equation (10), or through a shot 
noise intensity measurement represented by Pw,t and 
defined in equation (60). According to our semiclassical 


( 98 ) 


and quantum calculations, both are affected by photon 
bunching; therefore, these quantities must be correlated 
if our results are correct. Conversely, if we can establish a 
correlation between these quantities, the correlation may 
be used together with the well-known results for bunch¬ 
ing noise in direct photon counting to establish a lower 
bound for the bunching noise that must also be present 
in the shot noise measurements. In this section, I present 
a fully quantum treatment of these topics. 

We are interested in evaluating the correlation 

/ +00 

dvidv 2 dv 3 ( I{vi)l\v 2 )I{y3)) 

-OO 

rT/2 

xWi^W*^) dfie i27r(l/1 - y2)tl 
J-T/2 

( 101 ) 


fT/2 

x / dt 3 e i2 ™ 3t3 . 

J-T/2 























14 


We thus require the sixth-order moments of photon op¬ 
erators, 

{l(vi)I*(vi)I(v 3 )) = J dv{dv' 2 dv' 3 

x (b\v[)b(y[ + vi)tf(i/ 2 + v 2 )b(v' 2 )b ] (v 3 )b(v 3 + v 3 )) , 

to be evaluated as usual by computing the 3! = 6 op¬ 
erator pairings. As in section [7J the noise filters W{v\) 
and W*(v 2 ) allow us to ignore the DC terms in those 
variables; however, we must now retain DC terms for 1 / 3 . 
We can thus neglect permutations involving the cycles 
(1) and (2), which leaves only (12)(3), (123), and (132). 
We find: 


the variance of the photon counts Nt- Indeed, if X and 
Y are two random variables with zero mean and finite 
variance, the Cauchy-Schwarz inequality holds: 

(Ilf < (A 2 ) (y 2 ) , (103) 

which establishes a lower limit for the variance of X, 

( 104 ) 


Now set X = Pw.t — (CW.r) and Y = Nt — {Nt). From 
equations (12), (56) and ( |57| we have 

(Y 2 ) = T[T + S r (0)] • (105) 


(12)(3) = [f + S’r(z'i)] r$(iq - v 2 )S(v 3 ) 

(123) = J dv[ n{v[)n{v\ + v 3 + v[) 

X [1 + n(zq + v[)\ S(v 1 -v 2 + v 3 ) 
(132)= J dv[ n(v{) [1 + n(zq + v' x )\ 

x [1 + n(vi +v[ — v 2 )\ S(v 1 - v 2 + v 3 ) . 


Performing the integrations indicated in equation (101) 
gives 

7*+00 


/ -too 

dv 1 \W(v 1 )ff [f + 5 r (iq)] 

-oo 

/ +oo 

dv i \W(v x )\ 2 

-oo 

x Sr(^i) + J dv(n(v()n 2 (v{ + zq) 

/ +oo 

dv 1 \W(v 1 )\ 2 [f+ 5r(z/i) + 5 r (0) 

-OO 

+ / dv[n 2 (v{)n(v[ + v i) . 


The (12) (3) term just gives the product of averages 
(Pw.t) (Nt), so 

(■ P\v,tN t ) — (Pw.t) (Nt) 


C+oo 


= T 


dv i \W(vi)\ z {f + 5 r (0) + 25 r (z/i) 

J — OO 

+ J dv[n(v[) [n(v[) + n(v[ + v\)\ n(v( + 1 / 3 )| . 

Most of the terms vanish if we design the noise filter 
to reject the low-frequency noise as illustrated by the 
hatched region in Figure [31 the terms that survive are 

(P W . T N t ) — (P W ,t) (N t ) = 2BT [f + S'r(O)] , (102) 

using our standard definition of the shot noise measure¬ 
ment bandwidth B (equation [62]). We therefore see that 
the shot noise intensity Pw.t and photon counts Nt are 
indeed correlated, as we expect if both are affected by 
photon bunching. 


The value of the correlation given by equation (102) 
allows us to set a lower bound on the variance of the shot 
noise intensity Pw.t given the well-established results for 


Using the know n value of the correlation (XY) given by 
equation ( 102 ), we have 


<J 2 P = (X 2 )> (r P^-= 4B 2 T [f + 5 r (0)] . (106) 

Dividing by the square of the mean value (Pw.t) = 
2 BTT gives the fractional fluctuations 


4 s 1 , Sr(0) 

(Pw.t) 2 ~ fT f 2 T 


1 + 71 
fT 


(107) 


where we have used equations (73) and (74) in writing 
the second expression. Thus, using a fully quantum- 
mechanical calculation, we have demonstrated that the 
shot noise intensity measurement must suffer at least 
the same \J 1 + n sensitivity degradation due to photon 
bunching as does standard photon counting. Compari¬ 
son t o th e result of the semiclassical calculation, equa¬ 
tion ( [75| , shows that the correlation bound does not in¬ 
clude tne 1 / BT noise term associated with a finite band¬ 
width for the shot noise measurement. This is to be ex¬ 
pected: the finite-bandwidth noise does not influence the 
direct photon counts, represented by Nt, and is therefore 
absent in the correlation (Pw.tNt) ■ 

The extension to the case of two detectors is straight¬ 
forward but will be omitted. However, it is easy to see 
that the results above may be applied independently to 
each detector in the two-detector setup shown in Fig¬ 
ure [T] Thus, the shot noise intensity for each detector 
must be correlated with its photocurrent. Furthermore, 
the t wo photoc urre nts are correlated, as demonstrated by 
Hanbury Brown & Twiss (1956). Therefore, the (high- 
frequency) shot noise intensities of the two detectors 
must also be correlated, even though the shot noise itself 
is not: this distinction, between moments of t he f orm 
(l 2 I 2 ) vs. (/ 1 / 2 ), is elucidated further in section 10 


9. COMPARISON TO THE RESULTS OF LIEU ET AL. 

In this section, I comp are the res u lts of the previ¬ 
ous sections with those of Lieu et al. (2015) for both of 
the regimes they examine, corresponding to long sample 
times AvAt » 1 and short sample ti mes TAt < 1. Her e 
At is the single-sample time defined by Lieu et al. (2015); 
to avoid confusion, I use At instead of their chosen sym¬ 
bol, T, and instead reserve T = NAt to signify the total 
time duration of the measurement required for the ac¬ 
quisition and integration of N samples. The concept of 
sample time does not arise in my calculations since I as¬ 
sume continuous time integration; however, a connection 
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TABLE 1 

Comparison of symbols in this paper vs. those of|Lieu et al.| ([2015]). 


Quantity 

This paper 

Lieu et al. 

(as used here) 

Lieu et al. 
(original notation) 

Number of samples 

(continuous) 

N 

N 

Sample time 

(continuous) 

At 

T 

Total measurement time 

T 

T = NAt 

NT 

Optical bandwidth 

Ai/ 

Av = 1/r 

1/t 

Shot noise bandwidth 

B 

B At = 1/2A t 

1/2 T 

Photon arrival rate 

T = nAv 

noAv 

no/r 

Photon rate fluctuations 

S r (0) 

yTtUq Av 

Vnnl/T 


can readily be made since the sample time At defined by 
Lieu et al. sets the shot noise bandwidth B At = 1/2At 
associated with their measurement scheme. I make use of 
this correspondence to compare the two calculations for 
the same total measurement duration T, and find that 
the results agree in the long sample time regime but dis¬ 
agree for short sample times. Thus, my results directly 
contradict the claim of Lieu et al. that bunching noise 
may be avoided in the latter limit. It is important to note 
that for both regimes, the total measurement duration T 
can be chosen to satisfy A vT >> 1, as I have assumed 
for my calculations; indeed, long measurement durations 
are essential for astronomical observations since the sen¬ 
sitivity improves as l/vT. To aid in comparison of the 
results, and for ease of reference in the discussion below, 
the relevant quantities and symbols used to represent 
them in both papers are provided in Table [T] 

Lieu et al. (2015) present a fully quantum calculation 
of the shot noise fluctuations for the two-detector setup 
illustrated in Figure [I] They use a very similar quantum 
formalism for photon detection that differs only in minor 
and inconsequential detail. For example, their definition 
of the operator representing the detector output mea¬ 
sures photon power instead of photon counts as can be 
seen from their equation (5). Moments of photon oper¬ 
ators are calculated in the standard way, by combining 
operators pairwise, as is appropriate for thermal radia¬ 
tion. Lieu et al. consider the detector output averaged 
over some measurement time At corresponding to the 
quantity 

lAt = Xtj^ dtI{t) (108) 

in my notation. Lieu et al. focus on the difference of the 
outputs of the two detectors in Figure [l] 

4(1) =h{i)-h{t) , (109) 

and calculate both the second and fourth moments of 
1 

I A,At = J dtI A (t ) . (110) 

Their fundamental conclusions rely on evaluation of the 
mean and variance of the sum of N consecutive measure¬ 
ments of [/a, At] ) obtained over a total time duration of 
T = NAt. This quantity may be expressed as 

N 

Pn, lkd = an) 

k =1 


where 


^ pk At 

lA,At( k ) = TT / dtI A (t) 

J [k—X)At 


( 112 ) 


are the consecutive time-averaged samples of the pho¬ 
tocurrent difference I A (t). In contrast, I first apply an 
arbitrary linear filter to the photocurrent, 


AW) 


(t) = f 

J —( 


dt'W{t-t')I A (t') (113) 


and then study the mean and variance of the shot noise 
intensity integrated over time, 


Pa = 


f T/2 


dt 


-T/2 


C\t) 


(114) 


as illustrated in Figure [Sj 

Although the definitions of Pa,lkd and P A superfi¬ 
cially appear to be different, these two quantities are 
closely related, as illust rated in Fig u re [9] The averaging 
over At performed by Lieu et al. (20151 may be rep¬ 
resented by a particular (and inflexible) choice for the 
linear filter, namely a time window function: 

- {r' °js?. 

According to our definition (equation [62| , this filter has 
a bandwidth 


2B At — 


/ +oo 

-OO 

/ +oo 

-OO 

1 

At ’ 


du \W At { V f 
dt W At (t) 


we will also need 


/ +oo 

dv \W At {v)t = 

-OO 


2 

3At 


(115) 


(116) 


Note that this filter does not reject DC or low-frequency 
noise, but these are automatically rejected anyway by dif¬ 
ferencing the curre nts in a two - detec tor setup. Another 
distinction is that Lieu et al. (2015) perform the time 


integration operation as a discrete sum rather than a 
continuous integration: the output of the square-law de¬ 
tector is sampled at times tk = fcAf, and then summed, 
as represented by the dashed box in Figure[9] This choice 
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does not significantly affect the results, though the dis¬ 
crete sampling operation of Lieu et al. may result in a 
minor degradation in performance due to noise aliasing. 



F ig. 9. — Sig nal flow diagram for shot-noise detection as proposed 
by|Lieu et al.|(|2015). The ordering of the filtering and differencing 
operations may be interchanged. This scheme differs from that 
shown in Figure [8] in two ways: 1) the choice of filter is fixed and 
corresponds to boxcar integrator with time duration At; and 2) the 
(slow) time integration is not continuous but is instead performed 
in a discrete fashion using a sampler and summer (dashed box). 
The sampler operates at a rate 1/At and is synchronized to the 
filter. 


Equation (29) in Lieu et al. (2015) gives the second 
moment of one sample: 


/ 7-2 \ _ n o 

' A’At/ - ^ > 


(117) 


where we have omitted their factor of Wq so that the oper¬ 
ator represents photon flux. According to their equation 
(7), their symbol r is related to the optical coherence 
time and is inversely proportional to the optical band¬ 
width, r ~ 1/Av. Thus, translated to our notation, 


(^a,a t) ~ noAv2B At — 2TB At (118) 


where T = uqAv is the photon rate before the beam¬ 
splitter. Meanwhile, the corresponding equation for our 
observable (equation [88]) reads 

(P A ) = 2BTT . (119) 


These may be reconciled by using equation (114), 



lim 

r->o 


(Pa) 

T 


2TB . 


( 120 ) 


Thus, our results for the mean value of the shot noise 


intensity agree with Lieu et al. (2015) if we make the 
replacements r —> 1/Av and B -V B At = 1/2A t. 

We no w turn to the varia nce of a single output sam¬ 
ple of the Lieu et ah (20151 setup, which they calculate 
using a clever evaluation of the eighth-order moments 
of the photon operators in which most of the terms are 
discarded since they cancel in the two-detector scheme. 
Specifically, they calculate 


ct 1,LKD — (^A.At) — (^A,At) 


( 121 ) 


and express the result as a fractional variance in their 
equation (32) when the number of samples N = 1, 


Here F(x) is a smooth function that allows both the 
AvAt << 1 and AvAt >> 1 limits to be examined, 
and is derived under the assumption of a Gaussian spec¬ 
tral profile for the thermal radiation. Comparison of the 
single-sample variance with our results requires use of 
the latter limit because we assume AvT >> 1, where T 
is the duration of the measurement, for evaluation of the 
Fourier integrals. In this limit, F(x) —t and 


2 

a l,LKD 



2 + 


3^ 

AvAt 


1 

TAt ' 


(124) 


We may safely assume that consecutive samples are un¬ 
correlated in the AvAt » 1 limit, because the sample 
time At is long compared to the optical coherence time 
r. Thus, the f ractional variance for a sum of N samples 
(equation |111|) would be 


2 

a N, LKD 

77T 

yA ,At / 


1 

N 



AvAt 


1 

TAt. 


(125) 


Mea nwh ile, our result for the fractional variance (equa¬ 
tion 


93) reads 


'Pa 


(PaY 


1 

T 


2 f dv\W(v)\ 4 S'r(O) 


[Jdv\W{v)\ 2 


+ 


r 2 


2 f dvdv'\W{v)\ 2 \W{v')\ 2 S r (v - v') 
f 2 [fdv\W(v )| 2 ] 2 


(126) 


The first term is readily evaluated using equations (1151 
and ( 116|). The numerator in the third term must be 
evaluated in our chosen limit AvAt » 1, so M>'l « 
Av for the integrals; we obtain 


°Pa 

(PaY 


1 

T 


4 W(0) 

-At + 3 —Ta 


(127) 


For the Gaussian spectral profile used by |Lieu et al. 

polsl, 


t*+oo 


r = 


, , no 

av n(v) = - > tiqAv 


> — oo 

f*+CJO 


5r(0) = / dvn 2 (v ) = ^ n ° —y 

J -oo T 


nn^Av 


' —OO 

so our result reads 

Up 


(PaY 


4 At 
"3 ~T 


AvT TT 


(128) 


(129) 


Correspondence with Lieu et al. (2015) is obtained by 
letting the total measurement time T coincide with AT At, 
the time to obtain N samples, yielding 


2 

°1,LKD 


T 2 

2 A,A t 


t /At\ r 

= 2 + 3 —F — + -- 

At \ t J At no 


Translated into our notation, this reads 


'1. LKD 


T 2 

1 A,A t 


= 2 + 


3 F {AvAt) + 1 


AvAt 


TAt ' 


( 122 ) 


(123) 


(PaY 


N 


3y/n 


AvAt TAt 


(130) 


This expression reproduces t he th ree terms of the Lieu et 
al. result stated in equation (125), which is derived from 
their equation (32), apart from a somewhat smaller nu¬ 
merical factor on our first term which likely results from 
our use of a continuous integration over the measurement 
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time T instead of a sum of discrete samples taken every 
At as illustrated in Figure [9l Note that the second term 
of the Lieu et al. result confirms our S'r(O) term, which 
is the signature of photon bunching in the shot noise in¬ 
tensity. Thus, our results agree in the long sample time 
limit, AvAt » 1. 

In contrast, there is_a major disagreement in the short 
sample time regime, TAt < 1, which corresponds to the 
use of a wide bandwidth for measurement of the shot 
noise, Ba* > T = n qAv » Au as illustrated in Fig¬ 
ures [3] and [4] Lieu et al. fi nd F(x) ss \/x for x « 1, so 
their result (equation 1231 in this limit becomes 


2 

<7 l,LKD 


T 2 

J A,A t 


5 + 


1 


1 


TAt TAt 


(131) 


In other words, the goal is to choose At small enough that 
the Poisson term dominates the single-sample variance. 
Their principal claim, namely that bunching noise can 
be avoided, rests on the statement that a sum of N such 


samples, as given by equation (111) and illustrated in 
Figure [9j and acquired over a total measurement time 
T = NAt, would have a fractional variance 


Av, LKD 


1 & 


l.LKD 


1 


(-P/v,lkd) N 


t 2 

-A,At 


AT At 


= f ^ , (132) 


as would be expected if the samples were statistically in¬ 
dependent. Note that the limit A vT >> 1 as required 
for our calculation of the same two-detector scheme is 
reached simply by choosing N » 1/AuA t, so the com¬ 
parison's immediate through use of equations (93), (68) 
and (l75|): 


’ Pa 


1 1 + n 

{Pi) ~BT + TT 


(133) 


The second term dominates because we have assumed 
B > T. Thus, our calculation gives a sensitivity of (1 + 
n)/TT while Lieu et al. find 1/fT; our result includes 
the (1 + n) bunching penalty, while Lieu et al. claim it 
can be avoided. 

10. RESOLVING THE CONTRADICTION 

Section [9] shows that the detecti on sensit iviti es der ived 
in this paper agree with those of Lieu et al. (2015) for 
long sample times but disagree in the short sample time 
regime that is of central interest. In o btain ing their 
result for short sample times (equation 132), Lieu et 


al. assume statistical independence but do not actually 
prove this by computing the correlations between sam¬ 
ples, ([/ A At{k)f [Ia, At(0] 2 /- They instead state, af¬ 


ter their equation (32): “...data in non-overlapping time 
periods are uncorrelated, because the correlation func¬ 
tion (Hp(t)I!p( 0)) is proportional to a delta function...”. 
Translated to our notation, their statement relates to 

(IA, At (k) I a. At (l)) 

/>kAt AAt 

= TATA / dtl / dt2 (Ait^lA^)) 

{tAt) J(k-\)At J(l-l)At 

= hl Af (134) 


This demonstrates that lA,At(k) are uncorrelated as Lieu 
et al. claim, which is to be expected because the spec¬ 
tral density of /a(£) is white as illustrated in Figure |4j 
Nonetheless, this does not mean that the squares of these 
random variables, [lA,At(k )] 2 , are uncorrelated. 

Thus, we see that the fundamental claim of Lieu et 
al. is invalidated by a simple error, the fallacy of the 
converse. Suppose we have two zero-mean random vari¬ 
ables, X and Y. If they are independent, they must 
be uncorrelated, because {XY) = {X) (Y) = 0. How¬ 
ever, the converse is not necessarily true. If it were 
true, we could claim that when (XY) = 0, we must also 
have (A 2 F 2 ) = ( X 2 ) (T 2 ), which is the statement upon 
which the Lieu et al. result rests. A simple counterex¬ 
ample suffices: suppose the joint distribution of X and 
Y is given by 

f(x, y) = \ [5{x -y) + S(x + y )] g(x) , (135) 

where g(x) is a Gaussian distribution with zero mean 
and variance cr 2 . This distribution cannot be factorized 
into the form f(x,y ) = f x (x)f y (y), so clearly X and Y 
cannot be independent. We may readily compute 


(X) = J 

[ dxdy f(x,y)x = 

j dxg(x)x = 0 

(y) = J 

j dxdy f(x, y)y = 

[ dx^[x - x\g(x) =0 

{XY) = j 

j dxdy f(x,y)xy = 

\ [x 2 - x 2 ] g(x) = 0 

<*> = , 

j" dxdy f(x,y)x 2 = 

J dxg(x)x 2 = cr 2 

( y2 > = J 

[ dxdy f(x,y)y 2 = 

\ [z 2 + Z 2 ] g(x) = cr 2 

(X 2 Y 2 ) = 

[ dxdy f(x,y)x 2 y 2 

= 7 ; [x A + X 4 ] g(x) = 3cr 4 


Thus (X 2 Y 2 ) ± (A 2 ) (Y 2 ) even though (XY) = 0. 

It is quite easy to see that the samples l\ At (fc) must 
be correlated using a simple physic al arg ument. Consider 
the quantity defined in (equation (112): 


-1 rkAt 

lA,At(k) = ~ / 

1X1 J(k-l). 


dt [Ji(f) - I 2 (t)\ . 


(fe-l)Ai 

There are only three events that can occur with non- 
negligible probability when TAt << 1, corresponding to 
a short sample time: a) detector 1 receives a photon; 
b) detector 2 receives a photon; c) neither detector re¬ 
ceives a photon. These events correspond to values of 
lA,At.(k) = {+1/At, —1/At, 0}, respectively; therefore, 
I A a t(k) takes on the value of l/(At) 2 if either detector 
receives a photon, and zero otherwise. Thus, in the limit 
LA^ << 1 , the measu reme nt scheme proposed by Lieu| 
et al. (2015|) (equation |111| can be expressed as 


N 


P N ,LKD = ^ [lA,At(k )} 2 = —^ 
fc=l ^ ' 


(136) 


where Nt is the total number of photons received by 
both detectors over the course of a measurement of du- 
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ration T = NAt. Note that a single detector, replacing 
the two detectors and beamsplitter, would also have re¬ 
ceived Nt photons during such a measurement, so the 
statistics of the Lieu et al. observable P/v,lkd must 
be the same as those for N T , corresponding to photon 
counting with a single detector. The statistics of the lat¬ 
ter are well k now n to be affected by bunching, as stated 
in equation (13); thus the Lieu et al. claim that the 
bunching noise can be avoided is contradicted not only 
by the calculations presented in this paper, but also by 
the extensive experimental and theoretical work on pho¬ 
ton bunching over the past six decades. A more rigor¬ 
ous discussion is given in Appendix [EJ which provides 
a detailed quantum-mechanical calculation that demon¬ 
strates that the samples l\ At (fc) are indeed correlated, 


and that accounting for these correlations in the sen¬ 
sitivity calculation leads again to the standard photon 
bunching penalty, in agreement with the calculations for 
both one and two detectors presented in sections Em 
and [U 
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APPENDIX 

A. GENERATING FUNCTION FOR A POISSON PROCESS WITH A TIME-VARIABLE RATE 


Here we evaluate the generating function introduced in equation (45), 

oo fc M 

G N (s) = (e sNT ) = Yb £ (Vii-Vik) ■ 


k\ 

k =0 = l 


If the indices are all distinct, we may write 


(yh-yik) — (yii) ■■■ (Vik) — )...r(tj fc )Atj X ... 


(Al) 


because the y t are independent. When one or more indices repeat, we may use y™ = yi (for m> 1) to again obtain a 
product of distinct factors. We are thus faced with the problem of partitioning the set of indices {H ••■*&} into one or 
more groups, where the indices belonging to a group have the same value, and indices belonging to different groups 
have distinct values. The number of partiti ons of k ob ject s into p groups is given by the Stirling number of the second 
kind, S(k,p ), which are nonzero for p < k (Blasiak et al. 2007). We therefore write 


N 


Y Vii-Vik =Y S ( k ’ P '> Y 


Vii ■■■Vip 


(A2) 


*i—*fe=i 


p=0 


where the prime on the second sum indicates that the indices take on only distinct values. We may make the 
replacement 

Y Vii-Viv^p' Y 


Vii'-Vip • 


(A3) 
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by considering permutations of the indices. Taking the average, 


Y ( Vii-Vip ) 

i \^>^2 • • • ^ip 


ii ^ ^2 ■ • • ip 



rft p r(ti)...r(t p ) 



(A4) 


Use of the following identity for the Stirling numbers 


OO 

Y s ( k ’P)-]^y p = ex P iv ( eX -!)] 

k,p=0 


(A5) 


allows us to evaluate the generating function, 


Gn(s) 



exp [y (e s - 1)], 


(A6) 


where 


n = f dtT(t) . 

Jo 

B. DETAILED EVALUATION OF SHOT NOISE FLUCTUATIONS: SINGLE DETECTOR 


(A7) 


The required average in equation (661 may be performed by considering the partitions of the indices (see also 
Picinbono, Benjdaballah & Pouget 19707eqn. 2.16): 


(: ViViVkyi) = (yi) (%) (yk) (yi) (Ai) 

+ Sij (yi) (y k ) ( yi) + S ik ( yi) (yj) ( yi) + 5a ( yi) (yj) (y k ) (Bl, B2, B3) 

+ 5jk (yi) (yj) (yi) + Sji (yi) (yj) (y k ) + Ski fa) (yj) (y k ) (B4, B5, B6) 

+ 5ij5 k i fa) (y k ) + 6 ik 5ji fa) (yj) + 5a5 jk fa) (yj) (Cl, C2, C3) 

+ 5ij5 ik fa) fa) + SjjSu fa) (y k ) + 5 jk 5ji fa) (yj) (C4, C5, C6) 

+ 5 ik 5u fa) (yj) + 5ij5 ik 5ii fa) . (C7, Dl) 


We have neglected to subtract the correction terms such as ( yi) 2 (y k ) fa) because, as in section [ 4 J they contain an 
extra factor of At,; and therefore will vanish in the continuum limit. The number of terms of each parti t ion cl ass, 
here labeled A, B, C, and D, is (1,6, 7,1) and follows the sequence of Stirling numbers 5(4, k) (Blasia k et al.|200 7), as 
expected. Taking the continuum limit and evaluating the Fourier integrals gives the following terms: 
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I { yi ) i *[ y 2 ) Hy z ) i *{ v A ) 


r<5(i/i) + (5f(i/i) Y8(v 2 ) + 5 f*(z/ 2 ) r( 5 (v 3 ) + 5f (i/ 3 ) r<5(i/ 4 ) + < 5 f*(z/ 4 ) 


+ 

+ 

+ 


T 8 {vi — v 2 - 1/4) + 5 T{vi — i/2 — 1/4) fi 5 (z/ 3 ) + jf (i/ 3 ) 
fi 5 (i/i + 1/3 - 1/4) + ST (1/4 + 1/3 - 1/4) f( 5 (z/ 2 ) + < 5 f*(z/ 2 ) 

r<5(—1/ 2 + z/3 - i/4) + <5f*(i/ 2 -1/ 3 + i/4) r<5(i/i) + jf(1/4) 


(Al) 


+ 

Y8{yi - v 2 ) + 8T{yi - v 2 ) 

F(5(z/ 3 ) + tff(z/ 3 ) 

F5(z/ 4 ) + jf*(z/ 4 ) 

(Bl) 

+ 

r<5(z/i + z / 3 ) + <5f (z/i + z / 3 ) 

f 8{v 2 ) + (5f*(z/ 2 ) 

f(5(z/4) + 8Y*{vi 

)] (B2) 

+ 

T8(ui — z/ 4 ) + <5r(z/i — z/ 4 ) 

f5(z/ 2 ) + 5f*(z/ 2 ) 

F5(z/ 3 ) + <5r(z/ 3 ) 

(B3) 

+ 

T8{v 2 - z/ 3 ) + <5f (z / 2 - z/ 3 ) 

r^z/r) + <yf(z/!) 

f5(z/ 4 ) + 6 f* (z/ 4 ) 

(B4) 

+ 

T8{v 2 + z/ 4 ) + <5f (z / 2 + z/ 4 ) 

f(5(z/i) + <5f(z/i) 

f(5(z/ 3 ) + <5f (z/ 3 ) 

(B5) 

+ 

f(5(z / 3 - z/ 4 ) + <5r(z / 3 - z/ 4 ) 

F<5(z/i) + <5f(z/i) 

W(z/ 2 ) + 5r*(z/ 2 ) 

(B6) 

+ 

F<5(z/i - z/ 2 ) + ^r(z/i - z/ 2 ) 

f(5(z / 3 - z/ 4 ) + (5f(z / 3 - z/ 4 ) 


(Cl) 

+ 

r«5(z/i + z /3 ) + <5f (z/i + z / 3 ) 

F(5(z / 2 + z/ 4 ) + (5f*(z / 2 + z/ 4 ) 


(C2) 

+ 

T8(ui — z/ 4 ) + <5r(z/i — z/ 4 ) 

f<5(z / 2 - z/ 3 ) + <5f*(z / 2 - z/ 3 ) 


(C3) 

+ 

F<5(z/i — z/ 2 + z/ 3 ) + 8t(vi — 

z /2 + z/ 3 ) r<5(z/ 4 ) + <sr*(z/ 4 ) 


(C4) 


+ rj(i/! - i /2 + 1/3 - i/4) + (1/1 - i/2 + 1/3 - 1/4) • 


(C5) 

(C6) 

(C7) 

(Dl) 


Averaging over the stationary process <5r(f) now involves evaluation of its third-order and fourth-order moments. 
However, these higher-order mo ments are not fully specified by the second moment, which is determined by the power 
spectrum given in equation (571, because SY(t) is not guaranteed to be Gaussian. Indeed, that <5F(t) is not Gaussian is 
shown in Appendix [F| Nonetheless, <5r(f) may often be approximately Gaussian, and we proceed with this assumption 
recognizing that it may introduce small, detailed differences with the full quantum calculation. However, as could be 
anticipated, the term describing the sensitivity degradation due to photon bunching (labeled Clb below) involves only 
a second-order moment of <5r(t) and is therefore secure. 

For a Gaussian 8T(t), and omitting the DC terms, we find 


(1(^)1* (v 2 )I(v 3 )i*(v 4 )\ 

\ / y,<5r 

= S r (z/i)Sr(i/3)<>(z/i - v 2 )8(v 3 - 1/4) + Sr(z/i)S r (z/2)£(z/i - ^4)8^2 - 1/3) (Ala, Alb) 

+ S r (i/i)S r (i/ 4 )£(i/i + l/ 3)<5(*/2 + 1/4) + fS r (z/ 3 )<5(z/i - v 2 )5{v 3 - i/ 4 ) (Ale, Bl) 

+ tS T {y2)8{vi + v 3 )5(v 2 + i/ 4 ) + f5r(t/2)^(2/i - ^4)8^2 - z/3) (B2, B3) 

+ FSr(i/i)5(i/i - V4,)8 (v 2 - i/ 3 ) + f Sr(vi)8(vi + z/ 3 )6(z/ 2 + 1/4) (B4, B5) 

+ r5r(i/i)5(i/i — v 2 )8 {v 3 — i/ 4 ) + F 2 5(j/i — v 2 )8 {v 3 — i/ 4 ) (B6, Cla) 

+ r 2 5(isi + v 3 )8 {v 2 + 1/4) + r 2 J(i/i - va) 8 (v 2 - z/3) (C2a, C3a) 

+ [Sr(z/i - z/ 2 ) + S'r(z / i + z/3) + St{v 1 — z/ 4 ) + Sr (z/4) + Sr(z/3) (Clb - C5) 

+ Sr (z/2) + Sr(z/i) + f] 8{v 1 — v 2 + z/ 3 — z/ 4 ) (C6, C7, Dl) 


noting that the Al term gives three contributions, Ala - Ale due to the Gaussian pairwise evaluation of the fourth-order 
moment of <5F, while the factors in terms C1-C3 combine to give two contributions each, e.g. Cla and Clb. 

We now evaluate the second moment of the shot noise intensity measure in the limit of a long measurement time, 
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Ai/T >> 1. The result is 


( Pw,t) = j dvidv 2 dv 3 dv4: (^ 1 (^) 1 *{v 2 )i{y 3 )i*( 1 / 4 )^ ^ 


x W{v 1 )W*{v 2 )W{v 3 )W*{v 4 ) 


r+T/2 


dte i2n( ' Vl ~ l ' 2)t 


l-T/2 


r+T/2 


-T/2 


( ftV 27r(l ' 3 “ I ' 4)t ' 


_ rj-iZ 


dv\W (v)\ 2 Sr(y) 


2 T J dv\W {v)\ A S 2 (v) 


+ 2T 2 r / dv\W(y)\ 2 / dv'\W{y')\ 2 Sr{y') 


+ 4TT / ^|lT(^)| 4 S' r (^)+T 2 T 2 


+ 2TT 2 / d^|W r (i/)| 4 +TS , r(0) 


1 2 


e^|W»| 2 


^| 1 T ( i /)| 2 


+ 2T J dvdv'\W(v)\ 2 \W(v')\ 2 ST{v — v’) 
+ 4T J dv\W{v)\ 2 J dv'\W(v')\ 2 S T {v') 


+ TT 


dv\W(u )\ 5 


(Ala, Alb + Ale) 
(B1 + B6) 
(B2+B3+B4+B5, Cla) 
(C2a+C3a, Clb) 
(C2b+C3b) 
(C4b+C5b+C6b+C7b) 
(Dl) 


The terms proportional to T 2 sum to give ( Pw,t) ; the remaining terms proportional to T give the variance 


°"p ~ ( Pw,t ) — ( Pw,t) 

=t \2 / dv\W (v)\ A Sy(v) + 4T [ dv\W{v)\ A Sr(v) 


dis\W(is)\ 2 


+ 2f 2 / dz'|lT(i')| 4 + <Sr(0) 


+ 2 j dvdv'\W(v)\ 2 \W(v')\ 2 Sr{v — v') 
+ 4 J dv\W(v)\ 2 J dv'\W{v')\ 2 Sr{v') 


+r 


^| 1 T ( i /)| 2 


(Alb+c, B2+B3+B4+B5) 
(C2a+C3a, Clb) 
(C2b+C3b) 
(C4b+C5b+C6b+C7b) 

(Dl) 
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C. DETAILED EVALUATION OF SHOT NOISE FLUCTUATIONS: MULTIPLE DETECTORS 
Following the approach described in Appendix [Bj and omitting the DC terms, we have: 


ia 1 ( Z'l) la 2 C) I a 3 (■ V 3 )■Ki (■ V 4 ) 

= ^f al (l/ 1 )5f* 2 (l/ 2 )5f o3 (P3)5f* 4 (j/ 4 ) 

+ S a l,a2 

+ <^al,a3 
+ <^al,a4 
+ S a 2,a3 
+ S a 2,a4 
+ <5a3,a4 
+ <5al,a2l5a3,a4 
+ ^01,03^2,04 
+ <5ol,a4<5a2,a3 


5f a3 (u 3 )st:M 


TalS^! - V 2 ) + (Sr a i(l /1 - P 2 ) 
fal^l + V 3 ) + <5r a i(Pi + P 3 ) (5f* 2 (l/ 2 )<5f*4(p 4 ) 

r a i5(pi - P 4 ) + 5f al (i/! -1/ 4 ) <5f* 2 (^ 2 )c5f a3 (i/ 3 ) 

f a 2^(^2 - v 3 ) + <5f 02(^2 - P 3 ) <Sr a i(l/i)(5f* 4 ( 1 / 4 ) 

r a 2 <5(z/2 + ^ 4 ) + <5r a2 (p 2 +1/ 4 ) <5f a i(pi)<5f a3 (z/ 3 ) 

— P 4 ) + ^f a 3 (l / 3 — P 4 ) 

fai^i - ^ 2 ) + <5r a i(l/i - V 2 ) 

f'al'K^l + 1 / 3 ) + <5f a i(l/i + V 3 ) 

fal^l - V 4 ) + <5r al (l/! - P 4 ) 

+ <5al,a2<5al,a3 f a i<5(l/i — I / 2 + 1 / 3 ) + (JFal^l — V 2 + V 3 ) 

+ <Jol,a2^al,o4 fal^l - U 2 - V 4 ) + 5f a i(z/i - V 2 - U 4 ) ST a3 (v 3 ) 

4“ ^ol,a3^al,a4 |r a i<J(i / i + v 3 — v 4 ) + <5f al (p 1 + v 3 — v 4 ) I 8r* a2 (v 2 ) 

+ <5a2,a3^a2,o4 f a2 <J(-J/ 2 + P 3 - U 4 ) + ST* a2 (v 2 - V 3 + I/ 4 ) 5T ai^Vi) 

+ <5al,a2'5al,o3^ol,a4fal(5(l / l - V 2 + V 3 - V 4 ) + ST al (ui - V 2 + V 3 - 1 / 4 ) 


fo3^(^3 - Vi) + (Sr a3 (l/ 3 - 1 / 4 ) 

fo2^(f'2 + Vi) + jf * a 2 (v 2 + I/4) 
Ta2^(^2 - V 3 ) + ST* a2 (l2 2 - U 3 ) 
8Ki(Vi) 


(Al) 

(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 

(B6) 

(Cl) 

(C2) 

(C3) 

(C4) 

(C5) 

(C6) 

(C7) 

(Dl) 
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Inserting this result into equation (89) and evaluating some integrals, we find 


F abcd = T 2 J dv\W(u)\ 2 Cg\v) J dv'\W{v')\ 2 C^{v') 
+ T [ dv\W{v )| 4 [Cl + 


+ T 2 I dv\W{v )\ 2 / dv'\W(v')\ 2 [5 ab T a C^ d (v') + 5 cd t,-C^[y')\ 




+ Tj dv\W{v)\ A [ 5 ac T a C%\v) + S ad T a C£’(v) 
+6 bc T b C^(v) + 6 bd T b c£\v) 

+ T 2 S ab S cd f a f c J dv\W{v )| 2 
+ T [6 ac S bd r a r b + S ad S bc f a T c ] J dv\W{v)\ A 


'.(r)/ 


+ TS ab S cd Cg\ 0 ) 


dv\W{v)\ 2 


+ Tj dvdu'\W{v)\ 2 \W{v')\ 2 [8 ac 5 bd + 5 ad 5 bc ]C^ ] - v') 

+ T j dv\W{v)\ 2 j dv'\W{v')\ 2 [s ab 5 ac C ( f d \^) + S ab S ad Ci r c \F) 
+SacSa d C { J b \v') + S bc 6 bd C£\v') 


T TS ab S ac 6 ad T a 


1 2 


dv\W{v)\ 2 


The fluctuations of the noise intensity of the difference current are obtained by considering 


(-Pa) — Pnn + P 2222 + 4Pi2i2 + 2P 112 2 — 4F 1112 — 4F 22 i2 
= 2-Fiin + 2F1122 + 4T1212 — 8F1112 , 


(Ala) 
(Alb + Ale) 
(B1 + B6) 

(B2+B3+B4+B5) 

(Cla) 

(C2a+C3a) 

(Clb) 

(C2b+C3b) 

(C4b+C5b+C6b+C7b) 

(Dl) 


(Cl) 
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as can be seen by computing (equation [85| and making use of the symmetry of the 50/50 beamsplitter. The T 2 
terms are eliminated by subtracting the square of the mean, 


op A — (Ta) — {Pa ) 2 — 2-Fini + 2F1122 + 4 -F1212 — 8F1112 

= 4 T f dv\W{v)\ A [C ( ^\v)} 2 + STT! [ dv\W{v)\ A C { ^{v) 


+ 4TT 2 J dv\W(v)\ 4 + 2TC[ r 1 ) {0) 

+ 4 T J dvdv'\W(v)\ 2 \W(iy')\ 2 c[l\v - v') 

+ 8T J dvdv'\W{v)\ 2 \W{v')\ 2 C^{y') 

* 1 2 
dv\W{v)\ 2 


1 2 


dv\W{v)\ 2 


+ 2TT 1 


+ 4 Tj dv\W(v)\*[c£\v)] 


+ 2TC{P(0) 


2 

1 2 


dv\W{v){< 


+ 8T j dv\W{v)\ A [C ( ^\v)} 2 + 8Tt l j dv\W(v)\ 4 C[{\v) 

+ 4Tf\ J dv\W{v)\ A 

+ 4T J dudv'\W{v)\ 2 \W{v')\ 2 C { ^{v-v') 

-16 T J di/|W(i/)| 4 [C{P] 2 (i/) - 16Tfi J du\W{v)\ A C { ^{v) 

-8T J dvdv'\W{v)\ 2 \W{v')\ 2 C { ^ {v'). 


(2 F im ) 


(2 F, 


1122) 


(4 F, 


1212 ) 


( — 8F1112) 


The resulting sum is 


a 2 PA =+8Tf 2 / dv\W(v)\ A + 4TC[ 4 \o) 


dv\W{v)\ 2 


■8 T J dvdv'\W(v)\ 2 \W(v')\ 2 C { i X \v-v') 


2 TTi 


1 2 


dz/|W»| 2 


(C2) 


D. QUANTUM CALCULATION: EVALUATION OF EIGHTH-ORDER MOMENTS 
In this appendix, we evaluate eighth-order moments of the photon operators that are needed for a quantum- 


mechanical calculation of the sensitivity of a shot-noise measurement scheme, as outlined in section[7]and equation (95). 
As described in that section, there are 9 operator permutations that give nonvanishing contributions out of the 4!=24 
possibilities. We will not evaluate all of these terms but instead choose a few that are instructive, including the term 
that is responsible for the sensitivity degradation due to photon bunching. 

We start with the 2 2 permutations, (12)(34), (13)(24), and (14)(23): 


(12)(34) = (b\v[)b{v. 2 )) (b( u i + vi)tf{v' 2 + v 2 )) (& f K^K)) ( b ("3 + V 3 p(F 4 + 1 / 4 )) 
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Performing the indicated averages and integrations gives 

(12)(34) = J dv 4 dv' 2 dv' 3 dv' A {n(— i^)[fi(!/j + zq) + 1]<5 (i/( + iq — z^ — ^ 2 ) 
xn{v' 3 )8{v' 3 - v' 4 )[n(v 3 + v 3 ) + 1]<5(z^3 + v 3 - v' 4 - u 4 )} 

= J du' A dv 3 n(y' 4 )[n{y' 4 + zq) + l]5(zq - v 2 )n(v 3 )[n(v' 3 + u 3 ) + l]S(v 3 — u 4 ) 

= [f + 5 r (i2i)][f + 5'r(r' 3 )]<5(^i - v 2 )5(v 3 - v 4 ) . (Dl) 

Comparison with the semiclassical calculation detailed in Appendix |B] shows that we have reproduced the terms Ala, 
Bl, B6, and Cla; these become proportional to T 2 after the integrations over time and are related to the mean value 
of Pw.t rather than its fluctuations. Next, we evaluate 

(13)(24) = (b\u[)b{u' 3 + v 3 )) ( 1 b{u[ + zq)& f ( u' 3 )) {b\v' 2 + v 2 )b{v' A )) {b(v' 2 )b\u' A + v 4 )) . 

Upon averaging and integrating, 


(13)(24) = J dv[dv' 2 dv 3 dv A {n(v[)8(u[ — v 3 — v 3 )[n(v[ + zq) + 1]8{v' 3 — v' x — zq) 
xn{v' 4 )8(v' 4 — v' 2 — u 2 )[n(u A + v 4 ) + 1 ]8{y 2 — v A — zq)} 

= J dv' 4 dv' A {n(i/[)[n(i^ + zq) + l]5(zq + v 3 )n(v' 4 )[n(v' A + v 4 ) + 1 ]8{v 2 + v 4 )} 

= [f + S r (^i)][r + S' r (t' 2 )]<5(i / i + v 3 )5{v 2 + zq) . (D 2 ) 

Comparison with the semiclassical calculation shows that we have reproduced the terms Ale, B2, B5, and C2a; these 
become proportional to T after the integrations over time and are therefore related to the fluctuations of Pw,t- 
Similarly, 

(14)(23) = [f + <S r (zq)][r + S’ r (^ 2 )]<5(^i - v 4 )8(v 2 - v 3 ) (D3) 


corresponds to the semiclassical terms Alb, B3, B4, and C3a, which again are fluctuation terms since they are 
proportional to T. 

We now turn to the 4 1 terms: (1234), (1243), (1324), (1342), (1423), (1432). The terms that we have derived so 
far using the 2 2 permutations reproduce the semiclassical results of section pi and represent 1/BT noise, or terms 
that vanish if we choose a noise filter W(v) that rejects the bunching noise component at low frequency illustrated in 
Figure [ 3 J The 4 1 permutations are more interesting. We start with 


(1234) = (v[)b(v 2 )) (6 f (i4 + V 2 )b{v 3 + zq)) {b ] {v' 3 )b{v A )) (b(u[ + u 4 )b\v 4 + zq)) . 


Averaging and integrating, 


(1234) = J dv' l dv' 2 dv' 3 di'' 4 {n(z 2 })d(z/( — u' 2 )n(y 2 + v 2 )5{v 2 + is 2 — v 3 — zq) 
x n(v' 3 )8(v 3 - u'^[n{v' A + 1 / 4 ) + 1 \8{v[ + zq — z^ — zq)} 

= J dv[du 3 {n(v[)n(v[ + v 2 )n{v 3 )[n(v 3 + zq) + 1] 
x 8(u[ — v 3 + zq - v 3 )8{u[ — v 3 + v 4 ~ zq)} 

= J dv' 1 n{v' 4 )n(v[ + v 2 )n(v[ + zq — zq )[n(y[ + zq) + 1] 
x 8[y 1 - zq + zq - zq) . 
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This represents a contribution to the shot noise fluctuations given by 

r+o o 

( p w,t) { i234) = / dvidu 2 du 3 dv 4 W (v{)W* (v 2 )W (v 3 )W* (v 4 ) 


X -^(1234) (Fl, ^2, ^3, ^4)<j(^l - ^2 + ^3 - ^4) 


r +T/2 


dte i27l(vi ~ U2)t 


' — T/2 
/*+oo 


l‘+T/2 


I-T/2 


dt'e i2 


= T 


dvi dv 3 1 W(iq) | 2 1 W(u 3 ) | 2 

J —OO 

x J dv[n(u[)n(v[ + zq)n(z/( + zq — v 3 )[n(u[ + zq) + 1] 

This term did not appear in the semiclassical analysis, and presumably represents non-gaussianity of the photon arrival 
rate fluctuations which are expected in the quantum calculation as shown in Appendix IF] but are neglected in section [5] 
However, for this term to be appreciable, the noise frequencies zq and v 3 must be comparable to or smaller than the 
optical bandwidth Azq this term does not contribute if we choose a cutoff for W(y) that is well above A v as shown in 
Figure [3] Next is the pairing 

(1243) = (b ] {v[)b{ z4)) (h ] {v' 2 + v 2 )b{v' A )) (&K + v i)tf (i/ 3 )) (&(*4 + + v A )) ■ 

Averaging and integrating, 

(1243) = J dv^dv^dv'^dv^ {n(z^()(5(z^( — z4)n(z4 + v 2 )S(i , ' 4 — v' 2 — zq) 

x [n(v[ + z/i) + l]5(z4 - z^ - ^i)[rz(z / 3 + ^ 3 ) + I]<5(z4 + zq - v 3 - ^ 3 )} 

= J dv[n(i'' 1 )n(u , 1 + z/ 2 )[l + n(z/ 3 + zq)][l + n(u[ + zq + zq)] 
x 5(i/i - z/ 2 + zq - zq) . 

This term contributes a shot noise fluctuation given by 


(■ p w,t) 


p + OO 


(1243) 


= T 


dui dv 3 | W (zq) | 2 1 W (zq) | 5 


x J dz/ , 1 n(z/()n(z/( + zq)[l + n(u[ + zq)][l + n(v\ + zq + zq)] . 


Again, this term is small if we chose the high-pass filter cutoff frequency well above the optical bandwidth Azx Note 
that there is a contribution 


/• + OO 


dzqdzAj|VF(zq)| 2 |lF(zAj )| 2 / dzqrz(zq)n(zq + zq) 


r*+oo 


= T 


dv\\W (v \)\ 2 Sriyi) 


(*+00 


dv 3 \W (zx 3 ) | 2 


(D4) 


that reproduces the semiclassical term C7. We skip ahead and look at 

(1432) = (tf(v[)b(v[ 4)) (&K + + ^2)) (6(t' 2 )6 t (z/ 3 )) (6(z4 + iq)6 f (z/4 + zq)) 

Averaging and integrating, 

(1432) = J dv^dv^dv^dv^ {n(z^()5(z^( — z4)[l + n(u[ + v 4 )\5(y 4 + zq — v' 2 — v 2) 
x [1 + n(v' 2 )\S(u' 2 - v' 3 )[ 1 + n(v 3 + zq)]<5(z4 + zq - v 3 - zq)} 

= J dz^[n(z 2 ()[l + n(v\ + zx 3 )[1 + n(z^ + zq — zx 2 )][1 + n(v[ + zq)] 

X 5(zq — Zq + Zq — Zq) . 

The product expands to eight terms. The first term is 


dis[n(v[) 5 (vi — zq -)- zq — zq) = T< 5 (zq — v 2 + z / 3 - zq) 


(D5) 
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and reproduces the semiclassical term responsible for Poisson noise, Dl. Another term is 


J dv[n(v[)n(y[ +v\ — v 2 )8{v x - v 2 + ^3 — P4) = Sriyi — V2)8(vi — ^2 + ^3 — v 4) 


(D6) 


and reproduces the semiclassical term Clb in Appendix |B| 

For the case of multiple detectors, the photon operators are decorated with a subscript to indicate the detector that 
they correspond to. Thus, we are interested in 

(1432) = v 'i)) ( b a(v[ + vi)b\(v' 2 + u 2 )^ (bb(v2)bt{i/ 3 )) (b c {v 3 + vz)b\(v 4 + p 4 )) . 

Averaging and integrating, 


( 1432 ) = j dv^dv'zdv'zdv'i {B ad (v[)8(y' 1 - v' A )[8 ab + B ba (v[ + v x )\8{v[ + v A - v 2 - v 2 ) 
x [8 bc + B cb (v' 2 )\ 5 {v 2 - u' 3 )[8 cd + B dc (u' 3 + v 3 )\8(v' A + ^4 - p 3 - ^3)} 

= J dv'[B ad (i'' 1 )[ 5 ab + B ba (y' x + v{)\[8 bc + B cb (v x + v x — v 2 )][< 5 c d + B dc (i/[ + v A )\ 
x < 5 (^i - u 2 + V3 - v A ) . 

Again there are eight terms. The first term is 


S ab 5 bc 5cd J dv[B ad (v[)5(i' 1 - v 2 + v z - v A ) = 8 ab 8 ac 8 ad T a 8{v x -v 2 + v 3 ~ v A ) (D7) 

and reproduces the semiclassical Dl term in Appendix [C] that is responsible for Poisson noise. The interesting term is 

SabScd J dv[B ad {v[)B cb (v[ + u 1 - v 2 )5(vi -v 2 + v 3 - u A ) 

= S ab ScdC£?{vi - v 2 )8{v x - z/ 2 + p 3 - 1/4) (D8) 


and reproduces the semiclassical term Clb in Appendix [C] 

E. CORRELATION OF THE LIEU ET AL. SAMPLES: A QUANTUM CALCULATION 
This appendix presents a quantum-mechanical calculation that shows that the samples d AA 4 (fc) introduced in 

(2015) to be independent are in fact correlated; furthermore, these 
10 ton”bunching penalty. Let I a {t) be the output of detector a; here 


equations (111 112 1 and assumed by Lieu et al. 
correlations are shown to lead to the standard p 


a = 1 or 2. The corresponding quantum operator is 


I a {t) = J dv 1 dv[bl(v 1 )bM)e- i2 ^-^ t . 


We define the integral over the time interval [kAt, (k + 1) At] as 

r(fc+l)Ai 

I kAt 


[ { 

Ia,At(k) = j 
JkL 


dt I a (t) ■ 


(El) 


(E2) 


The Lieu et al. detection scheme involves summing the squares of the differences I%Ak) = I\,At{k) — I 2 ,At(k) of the 
two time-averaged and sampled outputs of a beanrsplitter-fed detector pair (Figures]!] and [9]): 


N—l 


S=J2 [lAt(k)] = At) 2 Pn, lkd 


k =0 


where Pjv.lkd is defined in equation ( 111 ). The mean value of a single sample is given by 


(E3) 


([4t(*0] 2 ) =G n (k)-G 12 (k)-G 21 (k) + G 22 (k) , 


(E4) 
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where 


G ab (k ) = (. I a ,At(k)Ib,At(k)) 

r{k+\)At 

= / dtidt2 

J kAt 

Ak+l)At 

= / dt±dt2 


J dv 1 di/ 1 dv 2 di/ 2 e- i 3 "te-^ tl e- i 2 "te-^ ta 


I kAt 


J dv 1 di/ 1 dv 2 di/ 2 e- av ^-^ tl e- i2 ^- ,/ '^ ta 


x {B aa (v- v'x)B bb (v2)5{v2 - v' 2 ) 

+ B ab (v{)5{v\ - v' 2 ) [<5 ab + B ba (u 2 )\ 5{v 2 - v[)} 


r a r;,(At) 2 + s ab T a At + 


dvB ab [y) 


(Aif ; 


(E5) 


we have made use of AvAt « 1 in approximating the third term. Note that the second term dominates in the 
short sample time regime, FAt << 1. Here = f dvB nn (y) is the photon rate for detector a; the quantity B ab {y) is 
introduced in section [3] through equations ( [27] ) , ( [36| , and |37| ). Therefore, 

(E6) 


([/i t (fc)] 2 }«(fi+f 2 )At = fAt. 


This is exactly what we expect given the discussion in section 10 for small At, the value of [d^ t (fc)] is unity if either 
detector receives a photon, and zero otherwise, and the probability of either receiving a photon is TAt. Thus, we 
conclude that mean value of the sum is 

(S) = NT At = FT (E7) 

where T = NAt is the total time duration of the measurement. 

To calculate the fluctuations of the sum S, we first define the quantity 

F ab cd(k ) = (/ a ,At(fc)^b,At(^)^c,At(0)Fd,At(0)) 


p(k+l)At f(k+l)At pAt pAt 

/ dti / dt 2 / dt 3 / dU(I a (ti)I b (t2)Ic(t 3 )I d (t4,)) . 

J kAt J kAt JO Jo 


We wish to evaluate the correlation 

CAt(k) = ([i d At (k)} 2 [i d At (0)Y 

We may easily express C/\ t {k) in terms of F abcd {k): 

CAt(k) = Tiin(fc) — F 2111 (k) — Ti 2 n(fc) + T 22 ii(A;) 

— F\\2i(k) + T 2 i 2 i(fc) + F 1221 (fc) — T 222 i(fc) 

— Fui2(k) + T 2 n 2 (fc) + F\2V2{k) — T 22 i 2 (fc) 

+ Fu22{k) — T 2 i 22 (fc) — Ti 222 (fc) + T 2222 (fc) . 

As usual, evaluation of F abcd {k) involves an eighth-order moment of photon operators, 


(E8) 


(E9) 


p(k+l)At p(k + l)At pAt pAt 

Fabcd{k) = dti dt 2 / dt 3 / dt 4 

J kAt J kAt Jo Jo 

poo 


X 


/0 


dv\ dv\ dv 2 dv' 2 d<v 3 d,i/ i dv 4 dv' i 

x ^^(vx-vTlt! e ~i 2 -K(w 1 -v! 1 )t 2 e -i2-n(v 3 -v' 3 )t 3 e —i2-K(v i -i' l i )t i 

Pairwise combination of the operators gives 4! = 24 terms. It is not difficult to show that the (12) (34) permutation 
gives 


FaiT 4 \k) = G ab (k)G cd ( 0) , 


(E10) 
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and inserting this result into equation (E9) gives 


C { At )m (k) = [Gn(0) - G 12 (0) - G 21 (0) + G 22 (0)]" 

= ([l d At(k )} 2 ) 2 ■ (Ell) 

This term will therefore contribute ( S ) 2 when calculating (S' 2 ), which will subtract out when we calculate the variance 
of S. As in Appendix [Dj the (12) (34) permutation contributes to the mean value rather than to the fluctuations. 

As before, the operator pairing corresponding to the (1432) permutation is responsible for the Poisson and bunching 
noise: 


A. 


(1432) 


abed 


/*(&+ 
(k) = / 

J kAt 


(k+l)At 


dti 


p(k+l)At 


pAt 


pAt 


dto 


dU 


dt± 


/ kAt 


x / dvidv' 1 dv2dv2di'3di'3di'idv' i 
Jo 

x e -i2n(i/ 1 -v / 1 )t 1 e ~i 2 Tr(i' 2 -v l 2 )t 2 e ~i2ir(v 3 -v' 3 )t 3 

X (ba(v[)bl{v 2 )) (bb{v' 2 )b\{y 3 )) lb c (v' 3 )b\{v 4 ) 


r(fc+l)Ai 


dt\ 


p(k+l)At 


pAt 


pAt 


dt 2 


dU 


dt^ / eh'idz/ 2 di' 3 <iz '4 


/ kAt 


/ ZcAt 


X e ~i 2 v(v 2 -i' 3 )t 2 

X B a d{ v l) [Sab + ^fjal^)] [<^bc + ^{,( 1 / 3 )] \8 c d + Bdciyi)] • (E12) 

This expression leads to eight terms; among these is the term that gives rise to the Clb contribution in the semiclassical 
and quantum calculations in Appendices |B] and [D] 


p(k+l)At /■ 4#c-t- 

F£T’ Clh (k) = dtj 

J kAt J kAt 


r (k+l)At 


r At 


r At 


dt 2 


dU 


dti 


dh'idv 2 dv^di'i 


x e -i2nv 1 (ti-t 4 ) e -i2Tvv 2 (t 2 -t 1 ) e ~i2m/ 3 (t 3 -t2) e -i2TTV4,(t4-t 3 ) 
x 8 a b8 c dB ac (vi )B ca (^3) ■ 

The integral gives S(t 2 — t\) while the V 4 integral gives S(t 4 — t 3 ); therefore 

r.(fe + l)At 


(E13) 


F. 


(1432),Clb 


abed 


p(k+l)At pAt p 
( k ) = dt 1 / dt 3 / 

JkAt. Jo Jo 


dvidv 3 


Note that 


x e -i 2 ^ ( t 1 -t 3 ) e -< 2 1 r« 8 (ts-ti)^ cdBac(l/l)Bco( ^ ) . 


B ac {vi)B ca {v 3 ) = -n(i/i)n(^ 3 ) 


(E14) 

(E15) 


regardless of the choice of indices. Furthermore, 

nOO 

/ dvidv 3 n(vi)n(v 3 )e 


— i 2 TTV\T g+i27TI/3T _ 


/»CXJ 

/ dvdu'n{v)n(v + v'y 2 ™' T 
Jo 

/»oo 

/ di/5r(i/)e iWT = A r (r) , 
Jo 


where Ap(r) is the Fourier transform of Sr{v) and represents the time autocorrelation function of the photon rate 
fluctuations. Note that Ap(r) decays on a timescale r ~ 1/Az/, and that A(0) = T 2 . Thus 

1 p(k+l)At pAt 

i'cf ),Clb (^) = -JabSed / dt 1 / dMr(ii - h) 

4 JfcAt JO 




r<5abl5cd(Af) 2 A r (fcAf) 


(E16) 


where the approximation holds because AtAv << 1 . The (1432) pairing also contributes a term that corresponds to 
Poisson noise, labeled D1 in the semiclassical calculation: 


F abfd ] '’ D1 (k) = ld ab S cd 5 ac S kfi Atf . 


(E17) 
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Of the sixteen terms in equation (E9), the only nonzero contributions for the Clb piece of the (1432) permutation 
come from Fnn, -F 1122 , -F 2211 , and -F 2222 , due to the S a bd c d factor; and all have the same sign. For the D1 piece, the 
additional S ac factor means that only Fim and F 2222 contribute. These two pieces give a contribution to C A t(k) given 
by 


c U432),Clb+Dl^ _ TAt§k o + (At) 2 A r (fcAt) . 


(E18) 


The second term in this expression shows that the quantities [l At ( k )] 2 are indeed correlated, in contradiction to the 
assumption of Lieu et al. (20151. The corresponding contribution to the variance of S (equation |E3|) is: 


4 = <S 2 > ^ <S > 2 

’ N-l 

Y c At (k~i) 


-<*5 ) 2 


k,l=0 
N-l 

= Y {TAt6 kl + (At) 2 A r [(k-l)At} } + ... 

k,l—0 


NT At. + N(Atj 2 


At 


+ ... 


where I made use of 


= TT{l+n) + ... 

1 /*+ OO 1 poo 

Y A r{( k - l ) At ] w At y A r (r)dT =-£j.J ri 2 (y)db 


At 


T = NAt, and T = nAv. Using equation (E7|, we find 


CT Af,LKD 


(Pn, 


lkd; 


(S)' 


1 +n 
TT 


(E19) 


(E20) 


(E21) 


This exp ressio n agre es w ith the othe r res ults presented in this paper (equation 133) but contradicts the fundamental 
result of Lieu et al. ( 2015| ) (equation 132). 

F. EQUIVALENCE OF QUANTUM AND SEMICLASSICAL APPROACHES 


In this Appendix, I use a straightforward extension of the arguments developed by Sudarshan 
the full quantum-mechanical calculation of photon bunching is equivalent to a semiciassicai ca. 


(1963) to show that 
dilation that makes 


use of a compound Poisson random process with a stochastically varying count rate. The equivalence is shown by 
comparing the generating functionals, defined as 


G[s] = ( exp 


p+oo 


dts(t)I(t) 


(FI) 


The semiciassicai and quantum-mechanical versions will be denoted by G ( . sc ' ) [s] and G^ qm ^[s], respectively. These 
generating functions fully encode the statistics of the photocurrent /(f); the statistics must be the same if G^ sc )[s] = 
G (qn b [si. 


A Poisson process with a deterministic time-varying rate F(t) obeys equation A6|: 

ptk + l 1 \ ( ptk+l 

'tk 


exp 


I(t)dt 


= exp 


T(t)dt 


(e s 1) 


If we make the time interval small enough, we may approximate 

rtk+1 


exp 


I{t)dt 


’t k 


exp {r(tfc)Atfc (e s — 1)} . 


J / y 


If the intervals [tk,tk+i\ span the region over which s(t) is nonzero, we may write 


/ +°0 _ ptk + 1 

dts(t)I(t) ~'Y / s(tk) / I(t)dt . 

-00 1, Jt k 


(F2) 


(F3) 


(F4) 
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Making use of the independence of the subinterval counts {y}, we have 


G[a 


exp 




r t k+1 

tk 
rtk+l 


I(t)dt 


=n( ex p s(t k ) f i(t)dt 
k \ L Jt k . „ 

~ exp |r(4) Atfc - l) } 

k 

= exp | ^2 r(4)A4 (e s(tfc) - l) 


and by taking the continuum limit we find 


G[s] = exp |y dt r(t) | . 


(F5) 


If we now allow the rate T(t) to be stochastic instead of deterministic, we must also perform an average over T(f). We 
obtain a formal expression for the semiclassical generating function by writing this average as a functional integral 


G (sc) [s] = J [dT(t))f[T(t)} exp | p dtr(i) (e'W - l) j , 


(F6) 


where /[r(f)] represents the probability density functional for the rate process T(f) (Ueda 1989) and [dr(f)] is the 
functional integration measure. 

We now show that the quantum generating function may also be written in this manner and obtain an expression 
for the resulting probabilitydensity /[r(i)]. The quantum-mechanical averages require traces over the thermal density 
matrix given by equation pi): 


G (qm) [s] = ^exp 
= TV 

= TV 

where s(u) is the Fourier transform of s(t) 


{ exp [/ 

Mi 


/■-too 

' dts(t)I(t ) 

— OO 

' /*+oo 


dt s(t)b\t)b(t) 


dv 1 dv 2 s{y 1 - v 2 )b\vi)b(v 2 ) 


/ +oo 

dts(t ) 

-OO 


i2-Kut 


(F7) 


(F8) 


In the following discussion, we will find it useful to switch between operators labeled by a continuous frequency index 
and a discrete approximation using 


pOO 

/ dvidv 2 s(vx - 1^2)^{pbp) -H- V' SijKbj 

Jo 


where 


rfi+Aui 


bi = 


\/A Vi , 


M) 


and similarly for b\, and therefore = Sij, while 


Sij = s(vi - Vj)y/ AviAvj 


(F9) 


(F10) 


(Fll) 


is a Hermitian matrix by virtue of s(—u) = s*(v). _ 

The coherent state representation is convenient for calculating G^ qm ^[s]. The coherent states (Glauber 1963) are 
given by 


I z) = exp |0) 


(F12) 
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and satisfy the normalization 


z') = exp ( 


Zi m = e~ 


(F13) 


where z represents the column vector with components {z^ and z' is its Hermitian conjugate, a row vector with 
components {z*}. The coherent states satisfy the overcompleteness relation 


where the integration measure is 


1 = J dp(z)e zfz | z) (z\ 

1 r d(ReZi)d(lmZi) 

Mz) = 11 - 


The thermal density matrix has a diagonal coherent-state representation 

p = det (./V^ 1 ) J dp{z)e~ z ^ z e~ z]N z \z) (, z\ 

where N is a diagonal matrix of mode occupation numbers with elements 

ddij — d f j — ~~ Sij . 

Thermal averages may be computed using this representation, 

Tr (Ap) = det (iV" 1 ) J dp{z)e~ zU ~e~ z]N ~ lz {z |A| z) 

The operator we are interested in has the form 

A = exp j ^2 b\ b j J = exp ( 'b^Sb ) 


(F14) 

(F15) 

(F16) 

(F17) 


(F18) 


where we use the vector notation for the photon operators in which b represents a column vector with elements {bi} 
and b 1 re presents its H ermitian conjugate. Coherent-state matrix elements may be evaluated using the normal ordering 
theorem (Blasiak et aI7||2007 1 

exp [b^Sb] =: exp [V (e s — l) b] : (F19) 

which gives a compact result for the quantum-mechanical generating function, 

G (qm )[s] = det (N- 1 ) J dp{z)e- zfz e- z ' N ~ lz {,z |: exp [tf (e s - l) b] :| z) 

= det (IV” 1 ) J dp(z)e~ z * z e~ z ' N z exp [V (e s — l) z\ {z\z) 

= det (N~ x ) J dp(z) exp [-«t (iV" 1 - e s + l) z] 
det (N- 1 ) 


det (iV _1 — e s + 1) 
1 


det [1 — N ( e s — 1)] 

= exp {—Trln [l — N (e s — l)] } . 

Here we have made use of the complex Gaussian integral ( |Negele fc Orland 1988]) 

exp [c^B~ 1 d] 


dp(z) exp [— z^Bz — c' z — z'd] = 


det B 


(F20) 


(F21) 


Equation (F20) resem bles o ther results for thermal radiation, e.g. those of Beenakker (1998). 

Although equation (F20) is a relatively simple expression for the quantum-mechanical g ener ating function, it is 
not easy to compare this result to our semiclassical generating function given by equation (F6). If we hold off the 
z-integration, we have 

G^ qm) [s] = det (TV" 1 ) J dp{z) exp [-z^N~ x z - z f (e s - 1 ) z\ . 
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In the continuum limit, the second term in the argument of the exponential is 

(e s — l) z = J dvdv'z*(v ) ( e s — l) vv/ z{v') . 


Now 


(e s - l) vy = s{v -v') + 7j\ J dv iK v ~ vi)K v i - v ') + ■■■ 

= /"■’('!) [^‘i] 

+ ^ J <ftidt2«(ti)s(i 2 ) J dv 1 e i2 ^- Ul)tl e i2 ^- l '' )t2 

= [ dtl S{tl) + ^yS 2 (tl) + ... 


0 i2n(i'— v)t\ 


= / dt 


e a(t) - 1 


i2n(iy— v')t 


Thus we obtain 


where 


z' (e s -\)z = j dtT(t\z) (e‘l‘>-l) 

T{t\z) = J dvdv’e i2ir< ' l '~'' s>t z* (y)z(v') . 


(F22) 

(F23) 

(F24) 


We thus conclude that the quantum-mechanical ge nera ting function may be written in a form identical to that of a 
compound Poisson process as expressed in equation 


(F25) 


G( qm )[s] = exp [—TrlnlV] J dfi{z) exp (—z^N 1 z) exp J dtT(t\z) (e sd> — 

= J [ dT(t)]f[T(t)\N} exp | J + °° dtT(t) (e s ^ - l) | 

= G (sc) [s] , 

provided that the probability density functional for the stochastic rate process is given by 

/[r(t)|AT] = exp [-TrlnlV] J dfi(z) exp (— z*N~ 1 z) <5 [r(£) — r(t|z)] (F26) 


and where r(t|z) is given by equation (F24|. Note that while z(y) has a Gaussian distribution with variance n(u), 
r(t|z) is a quadratic form of z{y) and therefore is not strictly Gaussian. 

It is not difficult to demonstrate that 


exp [—TV In AT] J d^i(z) exp (— z^N l z) z(vi)z*(v 2 ) = n(v\)5{v\ — u 2 ) 


(F27) 


while 


exp [—TrlnlV] J dfi(z) exp (— 2 ^ N x z) z(vi)z*(is 2 )z(v 2 )z*(i / 4 ) 

= n{v\)n{y 3 )8{v\ - v 2 )8{v 3 - u 4 ) + n(ui)n(v 3 )S(u 1 - v 4 )8(v 3 - u 2 ) . 
Thus, the mean of the equivalent stochastic rate process is 

(r(G)> = J drtrmnmmh) 

= exp [-TrlnlV] J dn(z)exp(-ziN~ 1 z)r(t 1 \z) 

= J dudv'^ v - v '^ tl n{y)8{v - v') 

= / dvn{y) = F , 


(F28) 
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which is the expected result. Meanwhile, the second moment is 

(T{h)T{t 2 )) = J du 1 du 2 du 3 du 4 e i2n ^ 1 - , ' 2 '> tl e i27r ^3-^)t2 


x [n(v 1 )6(i> 1 - v 2 )n(v 3 )8{v 3 - i/ 4 ) + n(i>i)6(i> 1 - v A )n{v 3 )5(y 3 - u 2 )\ 


= / dv\dv 3 


n{yi)n{v 3 ) + 


= r 2 + j dfSr [y) e i2nv( - tl ~ t2 ^ 


These results coincide with equations (56) and (|57|) . 


(F29) 





